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Introduction 

Fields  as  diverse  as  patient  safety  improvement,  body  armor  solution  design,  creation  of  new 
surgical  instruments  and  rehearsal  of  new  surgical  procedures  are  all  benefiting  from  the  creation 
of  computer  models  of  living  tissues.  These  models  have  applications  in  surgical  simulation 
systems  that  allow  new  doctors  to  experience  their  first  surgeries  without  risk  to  real  patients. 
They  can  be  implemented  in  finite  element  models  that  describe  the  interaction  between  an 
armor  solution  and  the  tissue  it  protects.  Doctors  and  engineers  can  test  new  instrument  and 
procedure  concepts  in  silico,  reducing  both  the  need  for  animal  testing  and  the  time  to  market  or 
clinical  application.  This  research  program  has  focused  on  creating  mathematical  models  which 
capture  the  behavior  of  tissues,  developing  instruments  to  measure  the  responses  of  tissues  to 
controlled  loading  and  deformation,  and  performing  the  experiments  to  determine  the 
characteristic  parameters  for  a  variety  of  non-load-bearing  organ  tissues  in  vivo.  Over  the  course 
of  the  initial  four  year  timeframe,  and  an  additional  extension  period,  we  have  progressed  in  all 
of  these  areas.  Our  models  have  evolved  from  simple  lumped-element  tissue  models  to  an 
advanced  seven  parameter  model  which  captures  the  complex  non-linear  and  poro-visco-elastic 
behavior  of  organs  such  as  liver,  spleen  and  brain,  and  we  have  created  mathematical  algorithms 
which  enable  real-time  interactivity  with  virtual  tissues.  We  have  applied  pre-existing 
instruments  and  developed  a  series  of  novel  systems  to  measure  the  responses  of  solid  organ 
tissues  to  small  and  large  deformations  and  over  a  wide  range  of  time  scales,  and  we  have 
created  an  organ  perfusion  system  which  enables  near  -in  vivo  experiments  to  be  performed  on 
the  lab  bench.  Finally,  by  applying  the  data  acquired  in  in  vivo  and  in  vitro  experiments,  we 
have  applied  mathematical  tools  such  as  inverse  finite  element  analysis  to  relate  the 
measurements  obtained  to  the  models  to  determine  sets  of  parameters  that  define  the  materials 
tested.  This  research  program  has  lead  to  a  variety  of  external  collaborations  and  provides  the 
basis  for  further  research  aimed  at  improving  mathematical  description  and  simulation  of  living 
tissues. 

Body 

The  following  sections  present  first  the  results  of  the  work  performed  during  the  no-cost 
extension  (NCE)  period  (“Year  5”),  followed  by  a  review  of  the  results  achieved  and 
contributions  made  over  the  course  of  the  whole  research  effort.  The  Year  5  results  are 
presented  in  the  context  of  the  tasks  described  in  the  NCE  request,  similar  to  the  direct 
comparison  with  the  Statement  of  Work  in  the  annual  reports. 

The  original  proposal  objectives  and  Statement  of  Work  are  presented,  and  the  review  of 
the  research  effort  is  presented  in  the  context  of  the  three  primary  areas  of  work:  testing 
instrument  and  system  development;  experimental  measurement  of  tissue  responses;  and 
mathematical  model  development  and  characteristic  material  parameter  determination. 

Year  5  Extension  Review 

The  NCE  request  was  made  to  make  use  of  funds  which  remained  at  the  end  of  the 
original  four  year  grant  period.  These  funds  enabled  us  to  continue  to  support  the  researchers 
and  in  particular  two  doctoral  candidates  for  nearly  an  additional  year.  The  elements  of  the 
research  program  that  we  proposed  were  the  following: 
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•  Perform  ex  vivo  testing,  per  the  protocols  developed  for  porcine  liver,  on  porcine  spleen 
and  kidney,  using  the  extracorporeal  perfusion  system  developed  earlier. 

•  Prepare  documentation  detailing  the  model  and  the  parameters  extracted  from  the  ex  vivo 
testing. 

•  Complete  development  of  a  balloon-catheter  device  intended  for  use  ultimately  on  human 
tissues,  and  perform  tests  with  it  initially  on  tissues  ex  vivo,  and  if  time  permits,  in  vivo 
on  porcine  models. 

•  Establish  sources  for  human  tissue  samples  and  prepare  protocols  for  performing  tests  on 
these  tissues. 

•  Complete  implementation  of  the  mathematical  models  developed  to  date  into  real-time 
forms  suitable  for  interactive  manipulation  of  virtual  organs. 

In  addition,  we  continued  work  on  elements  from  Year  4: 

•  Developing  algorithms  to  extract  3-D  strain  field  information  from  3-D  ultrasound  scans 
of  tissue,  applying  3-D  strain  data  to  improving  uniqueness  of  tissue  parameter 
determination. 

•  Improved  perfusate  “recipe”  (to  enable  in  vitro  testing  to  at  least  6.5  hours),  perfusion 
system  configuration,  and  motorized  indenter  hardware 

•  Maintained  and  extended  external  collaborations 

Ex  vivo  testing  on  additional  organ  systems 

Figure  1  shows  a  harvested  porcine  spleen  being  tested  with  both  the  large  deformation 
creep  indenter  and  the  motorized  indenter  system.  One  goal  of  this  research  was  to  begin  to 
create  an  atlas  of  tissue  properties  across  a  wide  range  of  tissue  types,  and  the  extension  of  the 
system  to  the  study  of  spleen  is  one  of  the  first  steps  on  this  path. 


Figure  1:  Creep  indenter  (left  rear)  and  motorized  indenter  (foreground)  testing  of  porcine  spleen. 

Typical  results  from  the  spleen  tests  are  shown  in  Figure  2.  In  Figure  2(a),  the 
indentation  strain  is  recorded  over  time  for  a  fixed  load,  exhibiting  a  similar  combination  of  short 
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and  long  term  responses  observed  in  our  porcine  liver  tests  (see  below).  Figure  2(b)  shows  the 
reaction  force  of  the  tissue  as  the  motorized  indenter  moved  at  constant  speed  to  an  indentation 
depth  of  approximately  9mm  (roughly  equivalent  to  30%  strain).  As  can  be  seen,  at  higher  rates, 
the  reaction  forces  due  to  poro-viscous  effects  are  larger  than  for  lower  rates.  In  addition, 
preliminary  results  show  that  the  spleen  is  lOx  softer  than  liver  tissue,  and  seems  to  be  slightly 
less  sensitive  to  the  presence  or  absence  of  perfusion,  however  this  result  remains  to  be 
confirmed  over  a  larger  range  of  experimental  conditions. 


■Force  vs  -Displ  for  ea  indentation:  spleenl  b 


Displacement  [mm] 


Figure  2:  (a)  typical  creep  response  of  porcine  spleen,  (b)  Series  of  ramp  responses  of  porcine  spleen. 


Documentation  preparation:  model  description  and  parameters  extracted 

Year  5  saw  the  completion  of  the  doctoral  thesis  and  successful  defense  thereof  by  Dr. 
Amy  E.  Kerdok,  as  well  as  the  printing  of  our  Journal  of  Biomechanics  paper  [Kerdok, 
Ottensmeyer  &  Howe,  2006].  The  pre-print  of  this  paper  was  included  in  the  2005  annual  report, 
and  Chapter  4  and  Appendix  B  of  the  thesis  are  exerpted  in  Appendix  A  of  this  document. 
Together,  these  two  documents  describe  the  models  developed  during  the  course  of  this  research. 
In  brief,  Figure  3(A)  shows  the  models  used  to  describe  the  linear  responses  from  testing  using 
the  TeMPeST  1-D  device  to  perform  small,  sinusoidal  deformations  over  a  range  of  frequencies. 
This  simple  model  was  sufficient  to  capture  the  visco-elastic  response  observed  and  demonstrate 
that  our  perfusion  system  (described  in  previous  annual  reports  and  in  the  overall  summary 
below)  created  conditions  much  closer  to  the  in  vivo  state  than  testing  tissues  in  an  unsupported 
in  vitro  state.  Figure  3(B)  shows  the  model  used  to  describe  the  long  duration  testing  performed 
with  our  creep-indentation  tester,  and  includes  elements  to  capture  both  the  short  time  constant 
behavior  and  the  slower  response  as  fluids  within  the  tissue  redistribute  over  time. 

The  models  of  Figure  3  are  useful  for  describing  the  observed  response  to  the  types  of 
indentation  performed  with  our  instruments.  However,  they  are  not  applicable  to  arbitrary 
loading  modes,  and  as  such,  a  constitutive  model  that  describes  the  behavior  of  liver  tissue  in 
general,  rather  than  in  these  particular  experiments,  is  necessary.  The  model  of  Figure  4, 
described  in  Appendix  A,  accounts  for  constant  volume  changes  in  the  shape  of  an  element  of 
tissue  (e.g.  shearing)  and  changes  in  its  volume  as  a  result  of  compressibility  and  porous  flow 
through  the  tissue.  Of  the  twelve  parameters  listed  in  the  figure,  only  seven  are  independent,  and 
the  ranges  of  values  determined  from  the  model  optimization  process  are  listed  in  Table  1.  The 
meaning  and  derivation  of  all  of  the  parameters  is  described  in  the  appendix. 

In  addition  to  the  full  non-linear,  poro-visco-elastic  model,  a  simplified  linearized  version 
is  also  described  in  Appendix  A,  in  the  excerpt  from  Appendix  B  of  Dr.  Kerdok’s  thesis. 
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Figure  3:  Lumped  parameter  models  of  liver  tissue,  (a)  First-order  visco-elastic  model  used  to  describe 
TeMPeST  1-D  measurement  responses,  (b)  Second-order  model  capturing  short  and  long  time  constant 
behavior  of  creep  response. 


l*0>  fourth  Al,  Gp,  Sp,  IV,  fw,  Y Q  Ks,  K(_jm,  Kq,  M 


Figure  4:  Physically-based,  non-linear  constitutive  soft  tissue  model,  (left)  Isochoric  (constant  volume) 
response  elements,  (right)  Volumetric  response.  Model  parameters  listed  for  each  portion  of  the  model, 
include  both  independent  and  dependent  variables. 


Table  1:  Constitutive  model  independent  parameter  ranges 


Ks 

[Pa] 

Gp 

[Pa] 

m 

sP 

[Pa] 

po 

[Pa] 

A,|_ 

K0 

[m4/Ns] 

Parameter  range 

6615 

40.3 

1.18 

2393 

1.35 

1.035 

0.86 

during  optimization 

25  048 

116.0 

2.50 

14  214 

28.96 

1.051 

11.2 

Balloon-catheter  minimally  invasive  testing  instrument  development 

An  ultimate  goal  of  this  research  has  been  to  move  towards  making  measurements  on 
human  tissues  in  vivo  as  a  final  validation  that  the  measurements  made  on  human  or  other  tissues 
in  vitro  with  our  perfusion  system  are  accurate.  With  the  exception  of  the  TeMPeST  1-D 
instrument,  which  was  originally  designed  for  minimally  invasive  surgical  use,  none  of  our 
instruments  are  suitable  for  intra-operative  use  in  humans.  To  address  access  to  human  tissues, 
we  observe  that  the  hepatic  vein  of  the  liver  has  an  extremely  thin  layer  of  epithelium  separating 
the  lumen  from  the  parenchyma  of  the  organ.  As  a  result,  endovascular  measurements  of  tissue 
response  would  allow  us  to  make  direct  measurements  on  the  parenchymal  tissue.  Further,  since 
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balloon  catheters  designed  to  occlude  blood  flow  through  the  hepatic  vein  to  measure  wedge 
pressure  already  exist  and  are  in  use,  we  investigated  the  possibility  of  using  a  balloon  catheter  to 
measure  the  response  of  parenchyma  by  inflating  a  fluid-filled  balloon  catheter  and  measuring 
the  pressure-volume  response.  This  measurement  is  analogous  to  the  force-displacement 
response  made  in  our  indentation  testing.  Further,  by  combining  these  tests  with  non-invasive  3- 
D  ultrasound  scans,  we  would  be  able  to  observe  the  geometric  boundary  conditions  and  the  3-D 
strain  field  of  the  tissue  during  testing. 

While  human  testing  has  not  been  approached  to  date,  we  developed  the  system  shown  in 
Figure  5  to  examine  the  feasibility  of  this  technique.  As  shown  in  Figure  6,  the  first  balloon 
expansion  shows  a  response  much  different  from  that  measured  in  air.  The  second  and  third  tests 
are  intermediate  between  the  two,  though  much  closer  to  the  free-air  test.  We  hypothesize  that 
this  may  be  due  to  tissue  damage  caused  during  the  first  test,  but  may  instead  be  a  normal 
response  of  the  balloon  to  initial  and  repeated  tests.  We  expect  to  investigate  the  use  of  a  device 
of  this  time  further  in  the  future.  Given  our  commitments  to  completing  the  documentation  and 
advancing  the  state  of  the  modeling  efforts,  in  vivo  tests  using  the  balloon  catheter  system  were 
not  pursued. 

Human  tissue  sourcing  and  protocol  preparation 

Despite  initial  intentions  to  progress  to  human  tissue  testing,  none  was  performed  during 
the  course  of  this  research  project.  In  future  research  programs,  we  will  begin  the  path  towards 
these  tests  earlier. 

At  the  present  time,  we  are  investigating  the  procedures  necessary  to  qualify  the 
TeMPeST  1-D  instrument  for  human  use,  and  determine  the  best  methods  of  sterilizing  the 
instrument.  This  investigation  was  commenced  in  particular  as  a  result  of  a  contract  to  construct 
a  duplicate  of  the  TeMPeST  instrument  for  Dr.  Yohan  Payan  at  Universite  Joseph  Fourier  in 
Grenoble,  France.  This  group  seeks  to  understand  the  biomechanics  of  the  soft  tissues  of  the 
human  oral  cavity,  and  its  effect  on  the  physics  of  speech  production. 

Real-time  implementation  of  soft  tissue  models 

A  parallel  effort  under  way  within  the  Simulation  Group  has  been  the  development  of  an 
open-source  software  architecture  that  will  enable  researchers  from  different  groups  to  combine  a 
variety  of  algorithms  for  real-time  soft  tissue  deformations,  collision  detection,  graphical  and 
haptic  rendering  and  other  elements,  into  the  same  simulation  system.  An  international 
consortium  has  been  formed  to  develop  this  framework  (called  SOFA:  Software  Open 
Framework  Architecture)  lead  by  Dr.  Stephane  Cotin,  with  members  at  INRIA  and  TIMC  in 
France,  ETH  Zurich  in  Switzerland,  the  Computer  Graphics  &  Visualization  Laboratory  in  South 
Korea,  and  Stanford,  Case  Western  Reserve  and  UC  Berkeley. 

One  part  of  this  program  evolved  from  the  early  work  on  the  Truth  Cube  (TCI)  described 
initially  in  the  Year  1  annual  report.  The  Truth  Cube  (described  in  the  Validation  section,  below) 
was  a  silicone  gel  cube  with  embedded  fiducial  markers  suitable  for  detection  during  CT 
scanning.  Large  deformations  were  applied  to  the  TCI,  and  the  motions  of  the  fiducials 
provided  information  about  the  real  internal  strains,  to  which  mathematical  models  of 
deformation  can  be  compared.  Research  conducted  by  the  SOFA  consortium  evolved  to  the 
point  where  more  complex  geometries  and  real-time  modeling  techniques  have  been 
implemented  and  are  being  compared  with  more  complex  real  object  standards  such  as  the 
cylindrical  beam  shown  in  Figure  7. 
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Figure  5:  (a)  Balloon  catheter  pressure-volume  tester  fluid  control  and  sensing  apparatus.  SM=servomotor; 
SM  pot=  servomotor  potentiometer;  Psensor=pressure  sensor,  (b)  Occlusion  balloon  catheter  inserted  into 
porcine  liver  hepatic  vein,  (c)  Fluid  filled  occlusion  balloon  catheter  in  hepatic  vein  in  sectioned  porcine  liver. 


Calibrated,  zerod,  Pressure-Volume  data 
compare  liver  (blue1  red,  grn)  vs  air  (blk) 

r  i  :  T  T" 


0  1  2  3  4  5  6  7 

Volume  [ml] 


Figure  6:  Pressure-volume  response  of  balloon  catheter  device  in  air  (pvcal2  for  calibration)  and  three  tests  in 
porcine  liver. 
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Figure  7:  Polymer  gel  beam  reconstructed  from  CT  scan  data,  and  simulated  beams  using  advanced  and 
linear  finite  element  modeling  techniques. 

The  non-linear,  poro-visco-elastic  model  described  in  Dr.  Kerdok’s  thesis  has  been 
provided  to  SOFA  project  personnel,  who  will  be  investigating  integration  of  a  simplified  and 
optimized  form  into  future  implementations  of  SOFA-compatible  deformation  modules. 

3-D  Optical  Flow  Field  Regularization  Using  Finite  Elements 

During  Year  5,  we  made  significant  progress  in  advancing  our  ability  to  utilize  the  data 
collected  from  the  3-D  ultrasound  system  that  was  integrated  into  our  motorized  indentation 
system  in  Year  4. 

In  the  field  of  computer  vision,  a  particularly  challenging  problem  is  determining  the 
motion  of  a  real  scene  from  a  sequence  of  still  images.  In  general,  the  problem  is  insoluble  based 
only  on  the  image  data,  and  require  a  technique  called  regularization  to  add  additional  constraints 
beyond  the  image  data  to  uniquely  determine  the  motions.  Mathematical  tools  exist  to  perform 
regularization,  but  for  real  deformable  objects,  a  physics-based  regularization  scheme  that  makes 
use  of  the  equations  describing  elastic  and  visco-elastic  solid-body  mechanics  has  been 
suggested  as  a  superior  technique. 

In  Year  5,  we  created  such  a  system  for  tracking  the  motion  of  tissue  in  three  dimensions 
using  the  3-D  ultrasound  scanner,  combining  algorithms  for  estimating  local  motion  with  a  finite 
element  representation  of  the  same  tissue,  thereby  improving  the  accuracy  of  the  calculated  strain 
field.  This  system  can  then  be  integrated  into  the  inverse  finite  element  modeling  algorithm  that 
we  developed  previously  to  improve  our  estimates  of  the  material  parameters  that  were 
summarized  above  from  Appendix  A. 

In  essence,  a  set  of  virtual  springs  is  attached  to  each  node  of  a  finite  element  model  of 
the  tissue.  The  ends  of  the  springs  are  displaced  away  from  the  node  by  an  amount  related  the 
magnitude  of  the  motion  at  that  location  in  the  optical  field.  The  stiffness  of  the  springs  is 
related  to  the  confidence  of  the  estimate  of  the  local  optical  flow;  where  image  quality  is  high, 
virtual  spring  stiffness  is  correspondingly  high,  and  the  finite  element  mesh  is  tightly  bound  to 
the  optical  (or  ultrasound  in  our  case)  image.  In  regions  of  low  image  quality,  the  virtual  springs 
are  looser,  and  the  finite  element  model  estimates  the  true  position  of  a  given  node  given  the 
combination  of  the  tissue  constitutive  model  and  the  boundary  conditions  imposed  by  the  other 
virtual  springs. 
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This  algorithm  is  illustrated  in  Figure  8,  which  begins  with  the  estimation  of  local  flow, 
and  ends  with  a  deformed  finite  element  model  from  which  can  be  calculated  the  motion  of  any 
point  within  the  tissue. 
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Figure  8:  Jordan  et  al.  algorithm  for  regularizing  optical  flow  calculations  through  the  use  of  finite  element 
analysis. 


Evaluation  of  this  technique  has  been  performed  on  synthetic  data  sets  created  from  static 
ultrasound  scans  which  were  deformed  using  a  forward  finite  element  model.  The  deformed 
ideal  model  was  degraded  with  differing  levels  of  Gaussian  noise  to  test  how  well  the  algorithm 
could  reconstruct  the  known  deformation  from  degraded  data.  An  example  of  the  sequence  of 
undeformed  ultrasound  scan  to  deformed  tissue  to  reconstructed  deformation  (uniaxial 
compression  in  this  case)  is  shown  in  Figure  9.  Partial  details  of  the  technique  are  described  in 
Appendix  B,  with  full  details  awaiting  acceptance  for  publication  (see  Refereed  Conference 
Papers,  Jordan  et  al.,  2007) 


Figure  9:  3-D  flow  field  reconstruction  using  finite-element  regularized  optical  flow  reconstruction  algorithm. 

Improved  perfusate  recipe  and  improvements  to  perfusion  system  and  motorized 
indentation  tester 

As  a  result  of  minor  observed  damage  in  histological  analysis  of  liver  tissue  as  a  result  of 
over-pressure  in  the  perfusion  system  and  the  original  “recipe”  for  the  perfusate,  we  continued  to 
improve  the  perfusion  system  in  Year  5.  Specific  changes  that  were  made  are: 

•  Install  a  double-boiler  heater  to  prevent  overheating  of  the  perfusate  supplied  to  the 
arterial  reservoir. 
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•  Reduction  of  the  volume  of  the  arterial  and  venous  reservoirs  to  minimize  exposure  to 
room  air. 

•  Alteration  of  arterial  and  venous  reservoir  elevations  to  reflect  porcine  blood  pressures 
rather  than  human  blood  pressures 

•  Alteration  of  the  perfusate  recipe  from  heparinized  Lactated  Ringers  Solution  to  a 
combination  of  five  liters  Dextrose  5%  Lactated  Ringers  Solution  to  one  liter  6% 
Hetastarch  (the  heparin  used  systemically  immediately  prior  to  sacrifice  and  during  the 
initial  flush  to  prevent  clotting  within  the  organ) 


Histological  analysis  of  tissues  maintained  under  the  new  perfusion  protocol  showed  less 
damage  even  after  6.3  hours  of  perfusion  time  (see  Figure  10). 
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Figure  10:  Histological  images  of  porcine  liver  with  H&E  staining  under  400x  magnification,  (a)  immediately 
after  sacrifice,  (b)  immediately  after  initial  flush  (note  lack  of  red  blood  cells  in  vessels)  (c)  after  6.3  hours  of 
perfusion  and  200  minutes  of  sinusoidal  indentation  testing.  Minimal  noticeable  difference  between  (b)  and 
(c). 


In  addition,  the  contact  head  of  the  motorized  indentation  tester  was  replaced  with  a 
system  that  adheres  to  the  tissue  sample  by  application  of  suction  (see  Figure  11)  through  a 
porous  polyethylene  filter  (to  avoid  deformation  of  the  tissue).  The  suction  tip  improved  the 
non-slip  contact  between  the  indenter  and  the  tissue  as  compared  to  the  original  version,  which 
relied  on  application  of  a  preload  force  to  maintain  contact. 
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Figure  11:  Motorized  indenter  suction  head  in  contact  with  porcine  liver. 

Full  research  program  summary 

Our  original  research  proposal  laid  out  a  program  of  essentially  three  main  components 
working  in  parallel:  create  tools  and  techniques  to  interrogate  the  mechanical  responses  of  soft 
tissues  to  controlled  forces  and/or  deformations  and  perform  experiments  in  vitro  and  in  vivo ; 
develop  mathematical  models  and  algorithms  to  capture  the  observed  tissue  responses,  determine 
the  values  of  the  characteristic  parameters  for  the  tissues  in  question  and  generate  optimized  real¬ 
time  implementations  of  the  models;  and  validate  our  instruments,  techniques,  and  models. 

The  following  sections  briefly  the  progression  through  the  five  year  effort  in  each  of 
these  areas  in  turn.  Additional  details  were  reported  in  the  earlier  annual  reports. 

Instrument  and  tissue  testing  protocol  development 

At  the  outset  of  this  work,  we  had  access  to  the  TeMPeST  1-D  instrument  that  was 
developed  by  Dr.  Mark  Ottensmeyer  through  research  previously  funded  by  the  Simulation 
Group.  Through  a  collaboration  with  Dr.  Daniel  Kalanovic  of  the  University  of  Tuebingen,  we 
performed  comparisons  between  indentation  and  rotary  shear  by  making  use  of  his  ROSA2 
instrument,  demonstrating  good  agreement  on  silicone  tissue  phantoms,  but  poorer  performance 
during  in  vivo  testing,  which  could  not  be  resolved  prior  to  Dr.  Kalanovic ’s  return  to  Germany 
(see  Figure  12). 

Dr.  Kerdok  and  the  visiting  scientist  Dr.  Takashi  Maeno  (Keio  University,  Japan) 
developed  the  T-needle  testing  system,  shown  in  Figure  13,  which  examined  the  response  of 
liver  parenchyma  deep  within  the  organ  by  compressing  tissue  between  two  cross  bars.  The 
complexity  of  interpreting  data  from  a  manually  driven  version  of  the  instrument  lead  to 
motorizing  it  to  improve  controllability,  however  even  with  this  improvement,  our  mathematical 
modeling  progress  was  not  yet  up  to  the  task  of  determining  parameters  from  the  non-uniform 
stress/strain  field  between  the  tips  of  the  T-needles. 

Simplifying  the  geometry  of  contact  with  the  tissue  lead  Dr.  Kerdok  to  construct  the 
creep  indenter  (Figure  14),  which  also  allowed  her  to  measure  the  long  duration  response  of 
tissue  to  a  constant  load. 
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Figure  12:  (left)  TeMPeST  1-D  minimally  invasive  tissue  compliance  testing  device,  (center)  ROSA2  torsional 
shear  tester  showing  rotary  shear  head  (surrounding  fixation  ring  removed),  (right)  ROSA2  performing  in 
vivo  measurements  on  porcine  liver. 


Vacuum 

source 


Barrel 

micrometer 

positioner 


■v. 


Figure  13:  T-needle  parenchymal  testing  instrument,  (left)  showing  T-needles  assembled,  (right)  showing  T- 
needles  embedded  within  porcine  liver  sample. 


Figure  14:  (left)  Creep  indenter  testing  porcine  liver  in  vivo,  (right)  Creep  indenter  and  TeMPeST  1-D  testing 
separate  lobes  of  same  perfused  porcine  liver. 

As  one  of  the  goals  was  not  merely  to  test  the  parenchyma  of  a  solid  organ,  but  also  its 
vessels  and  capsule,  Dr.  Anna  M.  Galea  applied  a  tactile  probe  originally  designed  to  detect 
lumps  within  breast  tissue  (Figure  15),  to  detecting  voids  within  liver  samples.  The  results  of  her 
modeling  efforts  (described  below),  showed  that  the  pressure  of  the  fluid  within  the  vessel  was 
more  important  than  the  mechanical  properties  of  the  vessel  walls.  In  addition,  Dr.  Kerdok 
performed  tensile  testing  experiments  directly  on  the  liver  capsule  after  carefully  harvesting  it 
from  the  fresh  organ  and  separating  it  from  the  underlying  parenchymal  tissue  (see  Figure  16). 

An  early  attempt  to  add  imaging  technology  to  our  data  acquisition  armamentarium,  and 
thereby  improve  our  measurements  of  tissue  strain,  was  the  development  of  a  video-based  knife- 
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edge  testing  system  (Figure  15),  that  could  extract  the  surface  contour  of  a  silicone  tissue 
phantom.  This  technique  was  superseded  by  the  next  series  of  measurement  systems. 


Figure  15:  (left)  Hand-held  tactile  probe.  Array  of  force  sensor  mounted  to  bottom  surface,  3-D  position 
tracking  device  mounted  to  top  of  handle,  (right)  Video  camera  view  of  knife-edge  tester  (tip  highlighted  in 
red)  indenting  silicone  gel  tissue  phantom. 


Figure  16:  Liver  capsule  testing  system,  (top  left)  Capsule  sample  holder,  (top  right)  Typical  force- 
displacement  response  showing  non-linear  response  up  to  ultimate  strength  and  failure  afterwards,  (lower 
left  and  right)  Images  of  porcine  liver  capsule  undergoing  tensile  testing. 

To  greatly  improve  our  ability  to  apply  controlled  loads  and  displacements  to  tissue,  Dr. 
Kerdok  lead  the  development  of  the  motorized  indenter  system  (based  on  the  TestBench  Series 
high  fidelity  linear  actuator,  Bose  Corp.,  EnduraTEC  Systems  Group,  Minnetonka,  MN),  which 
can  apply  larger  deformations  than  the  TeMPeST  instrument,  with  a  significant  fraction  of  the 
mechanical  bandwidth  of  that  device.  She  integrated  it  with  the  perfusion  system  that  was 
developed  in  parallel  with  the  creep  indenter,  and  later,  with  Petr  Jordan,  added  the  scan  head  of 
a  3-D  ultrasound  system  beneath  the  plate  on  which  the  tissue  samples  rest.  This  last  addition 
enabled  us  to  examine  the  3-D  strain  field  within  the  tissue  while  the  force  and  displacement  of 
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the  indenter  tip  were  simultaneously  being  recorded,  thereby  improving  our  knowledge  of  the 
relevant  boundary  conditions  (see  Figure  17). 


Figure  17:  (left)  Motorized  indentation  tester,  (center)  3-D  ultrasound  probe  mounted  underneath  testing 
platform,  (right)  Detail  of  ultrasound  probe  mounting. 

The  Normothermic,  Extracorporeal  Liver  Perfusion  (NELP)  system  is  a  very  significant 
advance  over  simpler  methods  of  preparing  tissues  for  mechanical  testing.  By  supplying  a  whole 
organ  with  a  suitable  solution  that  maintained  osmotic  pressure  and  osmolarity,  physiological 
temperature  and  crucially,  physiological  blood  pressure  by  elevating  reservoirs  of  physiological 
solution  to  provide  appropriate  hydrostatic  pressures,  we  were  able  to  demonstrate  near  -in  vivo 
test  conditions  in  the  in  vitro  environment  of  the  lab  bench  (Figure  18).  The  NELP  system 
enabled  extended  testing  (up  to  6.3  hours)  of  tissues,  a  period  significantly  longer  than  would 
have  been  possible  in  vivo,  and  testing  in  a  much  more  controlled  environment  that  can  be 
achieved  intra-operatively. 


u 

n 

n 

39  C 

. o 


Heater 


Pump 


Figure  18:  Normo-thermic  Extracorporeal  Liver  Perfusion  (NELP)  system. 
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While  our  research  focused  mainly  on  porcine  tissues  (liver  and  spleen  in  particular), 
through  collaborations  with  other  groups,  we  were  able  to  make  in  vitro  and/or  in  vivo 
measurements  on  porcine  muscle,  fat  (with  Dr.  Elisa  Konofagou,  currently  at  Columbia 
University)  and  brain  tissue  (with  Prof.  Simona  Socrate,  MIT),  murine  (rat)  liver  and  kidney 
(with  Cynthia  Bruyns  and  Dr.  Kevin  Montgomery,  Stanford  University),  rat  and  mouse  skin  and 
scar  tissue  (Dr.  Saja  Scherer  at  Brigham  and  Womens  Hospital,  Boston),  and  bovine  vocal 
tissues  (Dr.  Steven  Zeitels,  Massachusetts  General  Hospital,  Boston). 


Modeling,  parameter  identification  and  optimization  for  real-time  application 

Our  modeling  efforts  progressed  from  simple  one  dimensional  first  and  second  order 
linear  approximations  of  tissue  behavior,  through  application  of  linear  and  non-linear  finite 
element  modeling  techniques  to  account  for  tissue  testing  geometry,  to  a  physics-based 
constitutive  model. 

As  mentioned  above  in  the  summary  of  Year  5  document  preparation,  the  TeMPeST  1-D 
measurements  were  fit  to  a  first  order  linear  model,  and  the  creep  indenter  data  to  a  second  order 
linear  model.  The  parameters  extracted  from  these  models  were  used  to  demonstrate  the  very 
significant  differences  in  behavior  between  tissue  tested  in  vivo  and  that  tested  in  vitro  with  no 
special  treatment,  and  that  the  use  of  our  perfusion  system  restored  nearly  all  of  the  in  vivo 
biomechanical  characteristics  to  harvested  organs. 

The  tactile  probe  described  above,  had  its  data  analyzed  by  comparison  with  finite 
element  models  of  parenchymal  tissue  surrounding  a  pressurized  void  representing  a  blood- filled 
vessel  within  the  tissue  (see  Figure  19).  This  modeling  effort  demonstrated  that  the  tactile  probe 
could  be  used  to  identify  the  presence,  size  and  depth  of  large  vessels. 


Figure  19:  FE  model  of  tactile  probe  measuring  reaction  forces  of  tissue  with  presurized  void.  The  simulation 
is  repeated  with  the  tactile  head  moved  incrementally  across  the  surface  of  the  simulated  tissue. 

The  most  significant  advance  in  our  modeling  capabilities  resulted  from  a  collaboration 
with  Prof.  Simona  Socrate  (MIT),  who  had  previously  developed  a  model  for  describing  the 
properties  of  the  tissue  of  the  cervix.  The  model  incorporated  elements  to  describe  the  non-linear 
character  of  extra-cellular  matrix  (collagen,  elastin)  stiffness,  and  elements  to  describe  the  visco¬ 
elastic  behavior  of  the  cells  and  porous  flow  of  fluids  through  the  tissue  (see  Figure  4,  above). 
Dr.  Kerdok  adapted  the  model  for  use  on  our  liver  tissue  measurements,  and  created  the  inverse 
finite  element  modeling  techniques  that  enabled  us  to  identify  the  parameters  of  the  tissue’s 
constitutive  model  based  on  the  measurements  of  surface  displacement  and  applied  force.  Initial 
evaluation  of  the  model  was  done  using  data  collected  from  force-displacement  measurements 
made  on  biopsied  human  breast  tumor  tissue,  collected  previously  in  Dr.  Howe’s  laboratory. 
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Dr.  Kerdok  also  identified  shortcomings  of  the  inverse  technique,  in  that  with  the  limited 
force-displacement  data,  it  was  difficult  to  uniquely  determine  the  values  of  the  parameters.  Mr. 
Jordan’s  main  contribution  was  the  development  of  the  optical  flow  techniques  described  in  the 
Year  5  summary,  which  convert  noisy  3-D  ultrasound  data  into  high  quality  approximations  of 
the  true  tissue  motion  within  the  organ  (see  Figure  8  and  Figure  9).  By  applying  these  additional 
data  to  the  inverse  problem,  the  range  of  possible  characteristic  parameters  can  be  significantly 
narrowed,  and  this  effort  will  continue  on  after  the  research  efforts  described  here. 

One  of  the  main  applications  of  the  tissue  models  and  parameters  determined  in  this 
research  is  the  creation  of  real-time,  interactive  simulations  of  surgical  interventions  for  training 
or  pre-operative  rehearsal.  Dr.  Cotin  has  lead  the  effort  to  develop  a  software  architecture, 
originally  known  as  CAML  (common  anatomical  modeling  language)  when  conceived  in  late 
2000  and  early  2001,  now  in  the  form  of  the  Software  Open  Framework  Architecture  (SOFA). 
SOFA  is  a  software  architecture  which  can  allow  the  interaction  of  different  real-time 
representations  of  tissue  deformation,  collision,  visualization,  haptic  interaction,  and  so  on. 
Optimized  versions  of  our  non-linear,  physics-based,  poro-visco-elastic  model,  will  be 
implemented  in  the  SOFA  format  for  wider  distribution  and  application. 

Validation 

While  our  goal  is  to  create  accurate  models  of  the  deformation  of  materials  that  exhibit 
non-linear  stiffness,  anisotropy,  non-Newtonian  porous  flow,  visco-elasticity  and  potentially 
other  complex  phenomena,  we  must  also  verify  that  our  models  can  accurately  describe  simple, 
well  characterized  materials.  To  support  the  model  validation  effort,  and  provide  a  set  of  gold- 
standard  data  to  the  wider  deformable-object  modeling  community,  we  created  the  Truth  Cube 
and  its  successor  standard  test  objects. 

The  original  Truth  Cube  was  a  clear  silicone  cube  (GE  RTV6166),  8cm  on  a  side, 
constructed  in  layers  with  Teflon  beads  embedded  in  a  regular  grid  with  1cm  spacing  (see  Figure 
20).  The  Cube  was  subjected  to  loading  conditions  including  uniaxial  compression  and 
hemispherical  indentation,  while  being  passed  through  a  CT  machine  to  record  not  just  surface 
deformations,  but  also  the  locations  of  the  displaced  Teflon  beads.  By  segmenting  the 
undeformed  and  deformed  scans  of  the  Cube,  the  internal  strain  field  could  be  calculated.  As  a 
difficult  challenge  for  finite  element  and  other  modeling  techniques  is  accurately  calculating 
large  deformations,  we  applied  up  to  18%  uniaxial  strain,  and  30%  nominal  hemispherical 
indentation  to  the  Cube. 

The  resulting  data  sets  were  posted  on  a  dedicated  website 
(http://biorobotics.harvard.edu/truthcubel,  and  compared  with  finite  element  simulations,  using 
known  material  properties  of  the  silicone.  In  many  cases,  the  mathematical  models  did  not  even 
converge  to  a  solution  at  the  large  strain  levels. 
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Figure  20:  (left)  Mold  for  Truth  Cube  1  fabrication  and  completed  Truth  Cube  1.  (right)  CT  image  through 
midplane  of  Truth  Cube  1  undergoing  spherical  indentation. 


Subsequent  to  the  Truth  Cube  testing,  we  experimented  with  implanting  the  Teflon  beads 
into  a  whole  porcine  liver  and  CT  scanning  the  organ  in  a  relaxed  state  and  in  a  deformed  state 
typical  of  the  organ  retraction  that  would  be  performed  during  an  operation  like  a 
choly cystectomy.  This  proof-of-concept  experiment  showed  that  the  use  of  Teflon  fiducial 
markers  was  a  suitable  technique  for  examining  the  large  strain  behavior  of  whole  organs. 

Subsequent  to  the  testing  of  the  original  Truth  Cube,  collaborators  at  ETH  Zurich 
(Switzerland),  Nagoya  University  (Japan)  and  Brown  University  encouraged  the  development  of 
a  second  standard  object  that  could  be  shared  amongst  the  groups,  tested  under  their  own 
conditions,  and  compared  against  all  of  the  others  (the  original  Cube  was  destroyed  during  the 
final  spherical  indentation  test).  Dr.  Vincent  Luboz  (Simulation  Group)  lead  the  fabrication  of 
the  new  (somewhat  misnamed...)  Truth  Cube  2’s  with  Dr.  Ottensmeyer  (Figure  21).  Dr.  Marc 
Hollenstein  (then  a  Ph.D.  candidate  at  ETH  Zurich)  performed  a  detailed  analysis  of  the 
compiled  data,  providing  indications  of  the  scope  and  applicability  of  the  various  testing 
techniques,  and  of  different  methods  of  simulating  the  testing  methods  applied  (Hollenstein, 
2005). 

The  SOFA  consortium  became  aware  of  the  Truth  Cube  series  of  objects,  and  recently 
requested  the  construction  of  a  new  object  subject  to  large  deformations.  This  lead  Drs.  Luboz 
and  Ottensmeyer  to  construct  cylindrical  polymer  gel  beams  which  deform  under  their  own 
weight.  Early  simulations  shown  in  Figure  22  clearly  demonstrate  that  linear  models  (which  run 
much  faster  than  more  detailed  non-linear  models)  do  not  predict  the  deformations  of  linear 
materials  subject  to  large  strains  well,  even  under  static  conditions  (let  alone  conditions  in  which 
loads  and  displacements  are  time  varying).  This  study,  supported  by  separate  funding,  is 
ongoing. 
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Figure  21:  Truth  Cube  2  (TC2)  undergoing  torsional  resonator  testing  at  ETH  Zurich.  This  test  object  has 
randomly  distributed  Teflon  beads  embedded  within  it,  rather  than  the  regular  grid  of  beads  in  TCI.  It  is 
made  of  Ecoflex  0030  (Smooth-On,  Inc.,  Easton,  PA) 


Figure  22:  Comparison  of  Ecoflex  0030  beam  under  gravity  loading  (real  object  segmented  from  CT  scan) 
with  geometrically  non-linear  and  linear  FE  models. 

Key  Research  Accomplishments 

Each  of  the  annual  reports  listed  that  year’s  particular  accomplishments.  Those  of  Year  5 
are  listed  here  first,  followed  by  an  enumeration  of  the  most  significant  elements  of  each  of  the 
previous  years,  grouped  by  category. 
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Year  5  Key  Research  Accomplishments 

•  Extension  of  tissue  testing  protocol  to  porcine  spleen,  collection  of  first  sets  of  data. 
Early  estimates  show  spleen  to  be  1  Ox  softer  than  liver  under  similar  conditions 

•  Completion  of  analysis  of  liver  tissue,  including  refinement  of  non-linear,  poro-visco- 
elastic  physics-based  model  and  determination  of  ranges  of  characteristic  parameters 

•  Initiation  of  adaptation  of  model  to  realtime  implementation,  with  goal  of  integration 
within  a  tissue  deformation  module  compatible  with  SOFA 

•  Construction  of  balloon-catheter-based  endovascular  tissue  testing  device.  Calibration  of 
balloon  elasticity  and  initial  experiments  in  porcine  liver  hepatic  vein,  showing 
significant  differences  from  balloon  inflation  in  air. 

•  Development  of  solid-mechanics-based  regularization  technique  for  improving  the 
estimation  of  3-D  strain  field  within  tissue  from  3-D  ultrasound  scanning  data 

•  Modify  the  NELP  system  to  correct  the  hydrostatic  pressures  of  arterial  and  venous 
perfusate  supplies,  improve  the  controllability  of  perfusate  temperature,  add  a  suction 
head  to  ensure  non-slip  conditions  between  the  indenter  head  and  tissue,  and  improve  the 
recipe  of  the  perfusate  to  minimize  observable  cellular  damage  during  extended  perfusion 
experiments. 

Overall  Research  Program  Accomplishments 

Instrument  and  Testing  Systems 

•  Design  and  construction  of  T-needle  parenchymal  tester  and  motorization  of  the 
instrument  to  enable  precise  control  of  motion  trajectories 

•  Evaluation  of  rotary  shear  as  modality  for  soft  tissue  testing 

•  Design  and  construction  of  creep  indenter  and  motorized  indenter  for  large  deformation, 
long  duration  and  high  rate  testing  of  whole  organ  tissues 

•  Establishment  of  the  Normo-thermic  Extracorporeal  Liver  Perfusion  system  to  enable  in 
vitro  testing  of  tissues  with  ability  to  capture  near  in  vivo  mechanical  responses 

•  Integration  of  3-D  ultrasound  scanner  with  motorized  indentation  system  to  better 
characterize  the  3-D  strain  field  generated  during  indentation 

Modeling  and  parameter  determination 

•  Progression  from  use  of  linear,  lumped  parameter  models  to  describe  tissue  response  to 
fully  3-D  physics-based  modeling  techniques 

•  Examination  of  reaction  forces  from  tactile  probe  in  response  to  the  presence  of  fluid 
filled  vessels  within  soft  tissues  using  finite  element  modeling  techniques 

•  Adaptation  of  model  for  rubbery  polymers  and  cervical  tissue  for  application  to  perfused 
solid  organ  tissues 

•  Creation  of  inverse  finite  element  optimization  technique  to  determine  the  seven 
independent  model  parameters  which  characterize  a  given  tissue  within  the  model. 

•  Establishment  of  mechanics-based  regularization  scheme  to  accurately  determine  3-D 
strain  field  within  deformed  tissues 

•  Conception  of  common  representation  and  interface  specifications  for  collaborative 
development  of  medical  simulation  algorithms.  Initial  Common  Anatomical  Modeling 
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Language  (CAML)  concept  has  “spun-off’  into  the  international  Software  Open 
Framework  Architecture  (SOFA)  consortium. 

Validation 

•  Construction  and  testing  of  Truth  Cube  large  deformation  standard  object,  comparison  of 
finite  element  models  with  measured  deformations  from  CT  imagery 

•  Testing  of  Teflon  fiducial  technique  as  applied  to  a  whole  porcine  liver 

•  Collaborative  development  of  Truth  Cube  2  standard  objects,  and  comparison  of 
measurement  and  modeling  techniques  between  groups  in  the  USA,  Europe  and  Japan 

•  Construction  of  third  generation  standard  object  to  examine  large  deformation  of  linear 
material  beam  configuration  to  static  loading  due  to  gravity.  Initial  comparison  with 
current  state-of-art  real-time  models 

Reportable  Outcomes 

Manuscripts 

Journal  Papers 

Kerdok  AE,  Ottensmeyer  MP,  Howe  RD.  The  Effects  of  Perfusion  on  the  Viscoelastic 
Characteristics  of  Liver,  Journal  of  Biomechanics,  39,  2221-2231,  2006. 

Konofagou  EE,  Ottensmeyer  M,  Agabian  S,  Dawson  SL,  Hynynen  K.  Estimating 
localized  oscillatory  tissue  motion  for  the  assessment  of  the  underlying  mechanical  modulus. 
Ultrasonics,  42(1-9),  951-6,  2004. 

Kerdok  AE,  Cotin  SM,  Ottensmeyer  MP,  Galea  AM,  Howe  RD,  Dawson  SL.  Truth 
cube:  Establishing  physical  standards  for  real  time  soft  tissue  simulation,  Medical  Image 
Analysis,  vol.  7,  pp.  283-291,  2003. 

Ottensmeyer  MP.  TeMPeST  1-D:  an  instrument  for  measuring  solid  organ  soft  tissue 
properties.  Experimental  Techniques,  26(3),  48-50,  May/June  2002. 

Refereed  Conference  Papers 

Jordan  P,  Zickler  TE,  Socrate  S,  Howe  RD.  Mechanical  Regularization  of  Optical  Flow: 
General  Framework  Using  Finite-Elements.  Pending  submission  to  IEEE  Computer  Society 
Conference  on  Computer  Vision  and  Pattern  Recognition,  Minneapolis,  MN,  18-23  June  2007. 

Liu  Y,  Kerdok  AE,  Howe  RD.  A  nonlinear  finite  element  model  of  soft  tissue 
indentation,  in  S.  Cotin  and  D.N.  Metaxas,  eds.,  Proceedings  of  Medical  Simulation: 
International  Symposium  -  ISMS  2004,  Cambridge,  MA,  June  17-18,  2004,  Lecture  Notes  in 
Computer  Science  vol.  3078,  Springer- Verlag,  pp. 67-76. 

Ottensmeyer  MP,  Kerdok  AE,  Howe  RD,  Dawson  SL.  The  effects  of  testing  environment 
on  the  viscoelastic  properties  of  soft  tissues, ,  in  S.  Cotin  and  D.N.  Metaxas,  eds.,  Proceedings  of 
Medical  Simulation:  International  Symposium  -  ISMS  2004,  Cambridge,  MA,  June  17-18,  2004, 
Lecture  Notes  in  Computer  Science  vol.  3078,  Springer- Verlag,  pp.9-18. 

Galea  AM,  Howe  RD.  Liver  Vessel  Parameter  Estimation  from  Tactile  Imaging 
Information,  in  S.  Cotin  and  D.N.  Metaxas,  eds.,  Proceedings  of  Medical  Simulation: 
International  Symposium  -  ISMS  2004,  Cambridge,  MA,  June  17-18,  2004,  Lecture  Notes  in 
Computer  Science  vol.  3078,  Springer- Verlag,  pp.  59-66. 
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Bruyns  C,  Ottensmeyer  M.  Measurements  of  Soft-Tissue  Mechanical  Properties  to 
Support  Development  of  a  Physically  Based  Virtual  Animal  Model.  Dohi  T  and  Kikinis  R,  eds. 
Proceedings  of  MICCAI  2002,  Laboratory  Notes  in  Computer  Science  2488,  35-43,  Springer- 
Verlag,  Berlin.  25-28  Sept.  2002,  pp282-9. 

Kerdok  AE,  Cotin  SM,  Ottensmeyer  MP,  Galea  AM,  Howe  RD,  Dawson  SL.  Truth 
cube:  Establishing  physical  standards  for  real  time  soft  tissue  simulation  International  Workshop 
on  Deformable  Modeling  and  Soft  Tissue  Simulation,  Bonn,  Germany,  2001. 

Contributed  Conference  Papers  and  Abstracts 

Valtorta  D,  Hollenstein  M,  Nava  A,  Luboz  V,  Lu  M,  Choi  A,  Mazza  E,  Zheng  YP,  Cotin 
SM.  Mechanical  characterization  of  soft  tissues:  comparison  of  different  experimental 
techniques  on  synthetic  materials  Proceedings  of  the  4th  International  Conference  on  the 
Ultrasonic  Measurement  and  Imaging  of  Tissue  Elasticity,  Lake  Travis,  Austin,  Texas,  USA, 
October  2005. 

Jordan  P,  Kerdok  AE,  Socrate  S,  Howe  RD.  Breast  Tissue  Parameter  Identification  for  a 
Nonlinear  Constitutive  Model  in  Proceedings  of  the  BMES  conference  in  Baltimore,  MD,  2005. 

Kerdok  AE,  Howe  RD.  Characterizing  Large  Deformation  Behavior  of  Liver  for 
Surgical  Simulation,  in  Proceedings  of  BMES,  Baltimore,  MD,  2005. 

Kerdok  AE,  Howe  RD.  A  Physical  Basis  for  a  two  Time  Constant  Constitutive  Model 
for  Liver,  submitted  to  the  2005  ASME  Summer  Bioengineering  Conference,  June  22-26,  Vail 
Cascade  Resort  &  Spa,  Vail,  Colorado,  2005.  (honorable  mention,  PhD  student  poster 
competition) 

Kerdok  AE,  Jordan  P,  Liu  Y,  Wellman  P,  Socrate  S,  Howe  RD.  Identification  of 
Nonlinear  Constitutive  Law  Parameters  of  Breast  Tissue,  submitted  to  the  2005  ASME  Summer 
Bioengineering  Conference,  June  22-26,  Vail  Cascade  Resort  &  Spa,  Vail,  Colorado,  2005. 

Kerdok  AE,  Socrate  S,  Howe  RD.  Soft  tissue  modeling  and  mechanics,  28th  American 
Society  of  Biomechanics  Annual  Conference.  X-CD  Technologies  Inc.,  Portland,  OR.  poster 
235,  2004. 

Kerdok  AE,  Howe  RD.  A  Technique  for  Measuring  Mechanical  Properties  of  Perfused 
Solid  Organs,  ASME  Summer  Bioengineering  Conference,  Key  Biscayne,  FL,  2003. 

Kalanovic  D,  Ottensmeyer  MP,  Gross  J,  Buess  G,  Dawson  SL.  Independent  testing  of 
soft  tissue  visco-elasticity  using  indentation  and  rotary  shear  deformations.  Proceedings  of 
Medicine  Meets  Virtual  Reality  11,  Newport  Beach,  CA.  IOS  Press.  ppl37-143.  Jan  22-25  2003. 

Ottensmeyer  MP.  In  vivo  measurement  of  solid  organ  mechanical  tissue  properties. 
Society  for  Experimental  Mechanics  Annual  Meeting,  Milwaukee,  WI.  pp328-333.  10-12  June 
2002. 

Ottensmeyer  MP.  In  vivo  measurement  of  solid  organ  visco-elastic  properties. 
MMVR02/10.  Proceedings  of  Medicine  Meets  Virtual  Reality  02/10,  Newport  Beach,  CA.  IOS 
Press.  pp328-33.  23-26  Jan  2002. 

Kerdok  AE.  Soft  Tissue  Characterization:  Mechanical  Property  Determination  from 
Biopsies  to  Whole  Organs,  Whitaker  Foundation  Biomedical  Research  Conference.  August  9-12, 
2001,  La  Jolla,  CA.  (Poster) 

Cotin  S,  Dawson  S.  CAML:  a  general  framework  for  the  development  of  medical 
simulation  systems.  Proceedings  of  SPIE;  2000.  p.  294-300. 

Cotin  S,  Shaffer  D,  Meglan  D,  Ottensmeyer  M,  Berry  P,  and  Dawson  S,  CAML:  a 
general  framework  for  the  development  of  medical  simulation  systems.  Proceedings  of  the  SPIE 
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conference  on  Digitization  of  the  Battlespace  and  Battlefield  Biomedical  Technologies,  Orlando, 
2000. 


Related  Theses 

Amy  Elizabeth  Kerdok.  Characterizing  the  Nonlinear  Mechanical  Response  of  Liver  to 
Surgical  Manipulation,  Doctoral  Thesis,  Department  of  Engineering  and  Applied  Sciences, 
Harvard  University,  Cambridge,  MA,  2006 

Marc  Hollenstein.  Mechanical  Characterization  of  Soft  Materials,  Diploma  (doctoral) 
Thesis,  ETH  Zurich,  2005. 

Lectures,  Workshops,  Tutorials,  Posters 

Dawson  SL,  Ottensmeyer  MP,  Kerdok  AE,  Howe  RD.  Enabling  Technologies  for 
Advanced  Soft  Tissue  Modeling,  presented  at  PRMRP  Military  Health  Research  Forum,  San 
Juan,  PR,  30  April  -  4  May  2006.  Poster  and  oral  presentations. 

Ottensmeyer  MP,  Neumann  PF.  Physical  and  Virtual  Simulators  for  Medical  Training, 
Dartmouth  College,  Engineering  Science  05:  Healthcare  and  Biotechnology  in  the  21st  Century, 
23  May  2006 

Kerdok  AE,  Howe  RD,  Dawson  SL,  Socrate  S.  Mechanical  characterization  of  liver  for 
surgical  simulation.  Presented  at  Industrial  Outreach  Program.  Harvard  University,  2005.  (poster 
presentation) 

Kerdok  AE,  Howe  RD.  Characterizing  Large  Deformation  Behavior  of  Liver  for  Surgical 
Simulation,  HST  Forum,  Boston,  MA.  (poster  presentation) 

Jordan  P,  Howe  RD.  Identifying  the  Parameters  of  a  Nonlinear  Constitutive  Law  for  Soft 
Tissue  Using  Three-Dimensional  Ultrasound  Imaging.  HST  Forum,  Boston,  MA.  (poster 
presentation) 

Ottensmeyer,  MP.  Biomedical  Engineering  (and  stuff...),  Shad  Valley  Dalhousie, 
summer  science  program  for  high  school  seniors.  26  July  2004 

Dawson  SL,  Ottensmeyer  MP,  Kerdok  AE,  Howe  RD.  Enabling  Technologies  for 
Advanced  Soft  Tissue  Modeling,  presented  at  PRMRP  Military  Health  Research  Forum,  San 
Juan,  PR,  25-28  April  2004.  Poster  and  oral  presentations 

Kerdok  AE.  The  Effects  of  Testing  Environment  on  Soft  Tissue  Properties,  presented  at 
HST  Forum,  Boston,  MA,  2004. 

Ottensmeyer  MP.  Tools  for  measuring  soft  tissue  properties.  Workshop  on  Reality- 
Based  Modeling  of  Tissues  for  Simulation  and  Robot-Assisted  Surgery,  at  IEEE/RSJ  IROS 
2003,  31  Oct  2003 

Ottensmeyer  MP.  Measuring  properties  of  living  tissue  or  How  to  make  Virtual  Organs 
feel  right.  Dartmouth  College,  Engineering  013:  Biotechnology  and  Virtual  Medicine,  29  July 
2003 

Ottensmeyer  MP.  In  vivo  measurement  of  the  mechanical  properties  of  soft  tissues. 
McMaster  University,  Mechanical  Engineering  Seminar  Series,  16  June  2003 

Ottensmeyer  MP.  In  vivo  measurement  of  the  mechanical  properties  of  soft  tissues.  City 
College  of  New  York,  City  University  of  New  York,  Dept,  of  Biomedical  Engineering  Seminar 
Series,  15  Apr  2003 

Kerdok  AE,  Howe  RD.  Measuring  Parenchymal  Properties  of  Perfused  Solid  Organs, 
MIT  Health  Sciences  and  Technology  2003  Forum  Abstract 
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Ottensmeyer  MP,  Kalanovic  D,  Gross  J.  Comparison  of  indentation  and  rotary  shear  as 
modes  for  interrogating  soft  tissue  visco-elasticity,  SMIT2002,  Annual  Conference  of  the 
Society  for  Medical  Innovation  and  Technology,  Oslo,  Norway,  5-7  Sept  2002. 

Ottensmeyer  MP.  In  vivo  measurement  techniques  for  determining  mechanical  responses 
of  soft  tissues.  Tufts  University  Mechanical  Engineering  Seminar  Series,  10  Apr  2002. 

Ottensmeyer  MP.  In  vivo  measurement  of  solid  organ  stiffness.  Stanford  Workshop  on 
Surgical  Simulation,  Stanford  University,  Stanford,  CA,  20-22  June  2001. 

Ottensmeyer  MP.  Design  of  Devices  for  Tissue  Property  Measurement  and  Haptic 
Display.  Part  of  half-day  tutorial  on  medical  simulation  and  training:  “Simulating  minimally 
invasive  surgical  procedures  in  virtual  environments:  From  tissue  mechanics  to  simulation  and 
training”  Medicine  Meets  Virtual  Reality  2000,  Newport  Beach,  CA,  27-30  Jan  2000. 

Degrees  Completed 

Ph.D.  Thesis 
Kerdok,  A.E. 

Funding  applied  for  based  on  this  work 

•  Medical  Simulation  Trainin  Initiative.  PI:  Dawson  SL.  Sponsor:  RAD  II/CCC. 
10/1/2006-9/30/2011  (proposal  in  preparation) 

•  Collaborative  Development  of  an  Open  Framework  for  Medical  Simulation.  PI:  Cotin 
SM.  Sponsor:  Telemedicine  and  Advanced  Technology  Research  Center,  US  Army. 
1/1/2006-12/31/2006 

•  Evaluation  of  Wound  Healing  and  Scar  Formation  Using  Mechanical  Impedance  Test. 
PI:  Ottensmeyer  MP.  Sponsor:  CIMIT,  under  Department  of  the  Army,  cooperative 
agreement  no.  DAMD 17-02-2-0006.  Proposal  rejected. 

•  Construction  of  Duplicate  TeMPeST  1-D  for  Dr.  Yohan  Payan,  TIMC.  PI:  Ottensmeyer 
MP.  Sponsors:  Universite  Joseph  Fourier  &  Project  MIDAS,  Grenoble,  France. 
12/1/2005-2/1/2006 

•  Simulation  of  Behind-armor  Effects  of  Ballistic  Threats.  PI:  Ottensmeyer  MP.  Sponsor: 
CIMIT,  under  Department  of  the  Army,  cooperative  agreement  no.  DAMD17-02-2-0006. 
1/1/2004-  12/12/2004 

•  Development  of  Vocal  Tissue  Property  Measurement  Instrumentation.  PI:  Ottensmeyer 
MP.  Sponsor:  Massachusetts  Eye  and  Ear  Infirmary.  1 1/15/2003  -  1 1/15/2004 

Employment/Research  Opportunities 

Dr.  Amy  E.  Kerdok  is  now  performing  post-doctoral  research  with  Prof.  Simona  Socrate,  Dept, 
of  Mechanical  Engineering,  MIT.  This  is  an  extension  of  the  physics-based  modeling  research 
elements  of  this  program. 

Databases  Established 

Results  of  the  original  Truth  Cube  standard  object  CT  scans  have  been  archived  and  made 
publicly  available  at: 

http :  //biorobotics  .harvard,  edu/  truthcube 
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Collaborations  Initiated 

•  SOFA:  Simulation  Open  Framework  Architecture  (ongoing) 

USA:  Simulation  Group,  MGH;  Stanford  Bio-computation  Lab,  Stanford;  Case  Western 
Reserve  University,  UC  Berkeley.  Europe:  INRIA  (Alcove,  Evasion,  Epidaure  Groups); 
TIMC,  ETH  Zurich.  Asia:  Computer  Graphics  &  Visualization  Lab  (South  Korea) 

•  Simulation  of  Behind-armor  Effects  of  Ballistic  Threats  (ongoing) 

Simulation  Group,  MGH;  Institute  for  Soldier  Nanotechnologies,  MIT 

•  Truth  Cube  2  (ongoing) 

USA:  Simulation  Group,  MGH;  Brown  University.  Europe:  ETH  Zurich.  Asia:  Nagoya 
University. 

•  Validation  of  Harmonic  Motion  Imaging  of  Soft  Tissue  Lesions  with  Indentation  Testing 
(5/2003-2/2004) 

Simulation  Group,  MGH;  Dept  of  Radiology,  Brigham  and  Women’s  Hospital. 

•  Comparison  of  Rotary  Shear  with  Normal  Indentation  Soft  Tissue  Measurements 
(4/2002-10/2002) 

USA:  Simulation  Group,  MGH;  Germany:  Section  for  Minimally  Invasive  Surgery, 
University  of  Tuebingen. 

•  Microgravity  Rat  Dissection  Simulator  (12/2001  -8/2002) 

Simulation  Group,  MGH;  National  Biocomputation  Center,  Stanford  University;  Center 
for  Bioinformatics,  NASA  Ames. 

Personnel  supported  through  this  research  effort 

Researchers 

Massachusetts  General  Hospital: 

Dr.  Steven  L.  Dawson,  M.D.,  (PI)  Department  of  Radiology,  The  Simulation  Group, 

Dr.  Stephane  M.  Cotin,  Ph.D.,  The  Simulation  Group 
Dr.  Mark  P.  Ottensmeyer,  Ph.D.,  The  Simulation  Group 
Dr.  Xunlei  Wu,  Ph.D.,  The  Simulation  Group 

Harvard  University 

Dr.  Robert  D.  Howe,  Ph.D.,  (Co-I)  Department  of  Engineering  and  Applied  Sciences 

Dr.  Takashi  Maeno,  Ph.D.,  visiting  scientist,  currently  in  Department  of  Mechanical 

Engineering,  Keio  University,  Japan 

Graduate  students 

Department  of  Engineering  and  Applied  Sciences,  Harvard  University 
Dr.  Amy  E.  Kerdok,  Ph.D. 

Dr.  Anna  M.  Galea,  Ph.D. 

Petr  Jordan,  S.M. 

Conclusion 

When  this  research  program  was  proposed  in  2001,  the  field  of  measurement  and 
modeling  of  soft  tissues  in  vivo  was  considered  to  be  a  “black  hole”;  while  cadaver  and  animal 
tissues  have  been  measured  in  vitro  for  decades,  there  was  a  significant  absence  of  data  available 
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for  those  wishing  to  create  simulations  of  the  deformation  of  living  tissues.  As  discussed  in  the 
Introduction,  a  variety  of  applications  require  better  mathematical  descriptions  of  tissues, 
including  surgical  simulation  for  training,  procedure  planning  and  instrument  design  (revolving 
around  improving  patient  safety)  and  others  such  as  high  rate  impact  simulation  for  the  design  of 
body  armor  solutions  (with  the  expectation  of  reducing  the  incidence  of  trauma). 

This  program  specifically  addressed  the  parallel  problems  of:  creating  instrumentation 
and  protocols  to  perform  controlled  loading  and  deformation  experiments  on  tissues  in  vivo  and 
(properly  prepared)  tissues  in  vitro ;  developing  mathematical  models  to  interpret  the  data 
collected;  and  validating  the  simplified,  optimized  models  created  to  do  real-time  approximations 
of  the  fully  detailed  versions. 

We  examined  the  types  of  in  vivo  testing  devices  that  had  been  described  in  the  literature 
at  the  outset  of  the  program,  and  determined  that  none  met  the  needs  that  we  expected  in  terms  of 
providing  repeatability,  well  defined  geometry  and  range  of  tissue  strains  and  strain  rates 
achievable.  As  a  result,  we  initiated  a  program  to  create  novel  instruments  to  address  these 
issues,  progressing  from  the  existing  TeMPeST  1-D  instrument  to  the  motorized  indenter 
combined  with  a  3-D  ultrasound  scanner.  A  major  contribution  of  the  measurement  part  of  the 
research  was  the  establishment  of  the  perfusion  system  to  enable  extended  testing  of  whole 
organs  in  a  near  in  vivo  state.  The  perfusion  system  provides  a  more  well  controlled 
environment  than  the  abdominal  cavity  of  a  test  animal  in  the  operating  room,  and  enables  the 
use  of  organs  harvested  from  other  researchers’  approved  non-survival  experiments.  This  latter 
capability  reduces  the  need  to  sacrifice  additional  animals  only  for  the  purpose  of  harvesting  one 
organ. 

As  our  instrumentation  advanced,  our  mathematical  modeling  tools  evolved  as  well. 
Simple  linear,  first  and  second-order  models  sufficed  to  demonstrate  that  the  perfusion  system 
created  the  desired  environment,  and  that  unperfused  tissues  behaved  very  differently  from 
tissues  in  vivo.  These  models  are  not  constitutive  models  suitable  for  simulating  the  generic 
deformations  that  would  be  generated  in  a  useful  surgical  simulator,  and  as  a  result  we  adopted 
and  implemented  a  modified  version  of  a  physics-based,  non-linear,  poro-visco-elastic  model 
which  captures  many  of  the  observed  behaviors  of  tissues.  To  determine  the  values  of  the 
independent  parameters  that  characterize  a  given  organ,  we  implemented  the  model  as  part  of  an 
inverse  finite  element  optimization  algorithm.  This  algorithm  begins  with  trial  values  of  the 
parameters,  calculates  the  simulated  response  given  known  boundary  conditions,  and  produces 
an  output  which  is  compared  to  measurements  made  experimentally.  The  parameter  estimates 
are  updated  based  on  the  difference  between  real  and  simulated  data  until  the  two  are  in  close 
agreement.  The  use  of  this  algorithm  has  allowed  us  to  determine  a  range  of  possible  values  for 
the  seven  parameters  that  correspond  with  the  observed  real  responses. 

To  improve  the  utility  of  the  model,  we  wish  to  reduce  the  calculated  range  of  the 
parameters,  which  results  from  force  and  displacement  data  that  is  insufficient  to  yield  a  unique 
determination.  This  uniqueness  problem  has  been  addressed  through  the  development  of  the 
optical  flow  techniques  described  above,  which  are  now  ready  to  be  implemented  within  the 
inverse  FE  algorithm  to  narrow  the  range  of  the  characteristic  parameters. 

As  models  cannot  be  assumed  to  be  accurate  by  comparing  them  with  simulations  of 
themselves,  we  established  the  series  of  standard  test  objects  that  have  come  to  be  known  as 
Truth  Cubes.  These  polymer  gel  objects  with  internal  fiducials  have  been  scanned  under  a 
variety  of  loading  conditions,  and  serve  as  a  gold  standard  against  which  current  and  future  soft 
tissue  deformation  models  can  be  evaluated. 
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This  body  of  work  has  resulted  in: 

•  publication  of  nine  peer-reviewed  scientific  journal  and  conference  papers 

•  publication  of  thirteen  contributed  scientific  conference  papers 

•  two  doctoral  theses,  one  on  the  soft  tissue  modeling  and  measurement,  the  other  on  the 
Truth  Cube  series  of  objects  and  an  international  collaboration  to  compare  different 
measurement  and  modeling  techniques 

•  at  least  17  seminars,  poster  presentations,  workshop  lectures  and  other  research 
communications 

•  six  grant  proposals  or  contracts  stemming  from  the  expertise  established  in  instrument 
development,  soft  tissue  characterization  and  mathematical  modeling 

•  six  national  and  international  collaborative  efforts  studying  additional  tissue  types, 
different  deformation  domains  and  common  approaches  to  medical  modeling  and 
simulation. 

Addressing  the  perspective  of  national  needs,  and  in  particular  those  of  the  armed  forces’ 
medical  personnel,  this  basic-level  research  provides  a  portion  of  the  knowledge  necessary  to 
create  better  medical  and  surgical  training  technology  for  the  medic  through  to  the  senior 
physician.  More  advanced  simulation  technologies  are  being  developed  both  within  this  research 
group  through  other  R&D  efforts,  and  by  numerous  research  and  commercial  groups  around  the 
world.  Many  of  them  will  rely  on  the  kind  of  data  that  we  have  produced,  to  ensure  the  realism 
and  accuracy  of  their  systems. 

When  one  considers  that  the  VA  healthcare  system  treats  not  only  combat-related  injuries 
and  conditions,  but  the  normal  medical  needs  of  the  members  of  the  armed  forces  and  their 
families,  it  is  clear  that  any  useful  advance  in  basic  medical  technology  will  lead  to  better  care. 
As  quoted  in  the  original  proposal, 

“...military  medical  personnel  have  almost  no  chance  during  peacetime  to 
practice  battlefield  trauma  care  skills.  As  a  result,  physicians  both  within  and 
outside  the  Department  of  Defense  believe  that  military  medical  personnel  are  not 
prepared  to  provide  trauma  care  to  severely  injured  soldiers  in  wartime. ...” 

[Bauer  et  ah,  2000] 

The  kinds  of  training  simulators  enabled  by  this  research  can  address  the  need  to  properly 
prepare  medical  personnel  for  the  conditions  that  they  will  face  both  in  the  field  and  in  the 
operating  room. 

For  the  future,  we  expect  to  continue  this  research  with  follow-on  funding  from  sources 
including  the  ones  listed  above.  Immediate  applications  of  the  data  collected  so  far  include: 
creation  of  improved  polymer-based  materials  to  create  physical  simulators  that  feel  more 
realistic  to  the  trainees;  implementation  of  the  models  into  the  SOFA  system  so  that  the  national 
and  international  partners  can  improve  their  models  of  various  organ  systems;  and  disseminating 
a  “recipe  book”  detailing  the  process  for  adapting  the  technology  and  models  developed  here  to 
researchers  studying  other  tissues  and  organ  systems. 
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Chapter  4 

Nonlinear  Physically-Based  Constitutive  Model  and  Parameter 
Identification 

The  mechanics  of  a  material  are  governed  by  the  behavior  of  its  constituents.  A 
constitutive  relation  determines  the  state  of  stress  at  any  point  in  a  body  in  response  to 
any  arbitrary  history  of  motion  (Holzapfel  2000).  Knowing  the  structure  and  the  constitutive 
relationships  that  describe  the  material  behavior  of  a  composite  body  (i.e.  an  intact 
organ),  the  response  of  the  body  to  any  prescribed  loading  condition  (i.e.  surgical 
manipulations)  can  be  obtained.  The  challenge  is  to  define  constitutive  relations  that  are 
both  simple  enough  to  be  easily  implemented  and  realistic  enough  to  accurately  describe 
the  material’s  behavior  across  different  loading  conditions.  Currently,  the  lack  of  elegant 
and  robust  constitutive  relations  for  soft  biological  tissues  hinders  the  development  of 
predictive  biomechanics  necessary  for  medical  simulation  (Fung  1993;  Humphrey  2002). 

Humphrey  proposes  a  five-step  procedure  for  developing  a  constitutive  relation 
(Humphrey  2002): 

1 .  Determine  the  general  characteristics  of  interest 

2.  Establish  a  theoretical  framework 

3.  Identify  the  specific  functional  forms  of  the  constitutive  relations 

4.  Calculate  the  values  of  the  material  parameters 

5.  Evaluate  the  predictive  capability  of  the  final  relation 

The  results  and  analysis  of  the  data  presented  in  Chapters  2  and  3  summarize  the 
general  mechanical  behavior  of  the  liver.  The  nonlinear,  viscoelastic,  strain-rate 
dependent  response  to  large  deformations  was  found  to  be  sensitive  to  its  geometric  and 
physiologic  boundary  conditions.  These  results  are  used  to  identify  the  functional  form  of 
the  constitutive  relations  under  the  theoretical  framework  of  continuum  mechanics.  The 
finite  element  method  is  used  to  solve  the  inverse  boundary  value  problem.  An  iterative 
approach  minimizes  a  prescribed  objective  function  to  identify  the  material  parameters  of 
the  model  from  the  multiple  ramp  indentation  and  stress  relaxation  data.  Finally,  the 
predictive  capability  of  the  model  and  its  parameters  are  validated  using  stress  relaxation, 
multiple  ramp  indentation,  and  the  creep  indentation  data. 
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4.1  Physiologically  Derived  Constitutive  Model 

In  Chapter  3  the  liver  was  described  as  being  mainly  comprised  of  a  tough 
collagenous  capsule,  vast  cellular  parenchyma,  and  pressurized  fluid-filled  vessels.  In  this 
chapter,  the  contribution  of  each  of  these  constituents  is  accounted  for  while  developing  a 
model  that  allows  for  the  cooperative  nature  of  the  tissue  response. 

The  multi-physics  based  continuum  model  developed  here  relies  on  the  previous 
modeling  work  of  other  soft  tissue  researchers.  Socrate  and  Boyce  developed  a  similar 
model  for  the  large  strain  behavior  of  articular  cartilage  (2001).  Febvey  further  developed 
their  model  to  describe  the  mechanics  of  cervical  stroma  (2003).  Both  of  these  solid  soft 
tissue  structures  were  modeled  as  hyperelastic  collagen  and  viscoelastic  proteoglycan 
networks  in  parallel  with  a  pressurized  interstitial  fluid  network.  Using  a  similar 
approach,  the  liver  is  modeled  as  a  homogenous,  initially  isotropic,  fluid-filled  structure, 
containing  a  hyperelastic  collagen  network  in  series  with  a  viscoelastic  cellular  network. 
A  third  porous  network  represents  the  volumetric  response  due  to  the  fluid  flow.  Finally, 
the  tough  outer  capsule  structure  is  modeled  as  a  hyperelastic  collagenous  membrane 
similar  to  Prevosf  s  treatment  of  the  chorioamnion  (fetal  membrane  sac)  (2006). 

4.1.1  Motivation 

The  extracellular  matrix  (ECM)  of  soft  tissues  regulates  the  cell’s  shape, 
orientation,  movement,  and  overall  function  (Park  and  Lakes  1992).  ECM  is  comprised  of 
elastin  and  collagen  proteins  in  a  hydrated  ground  substance.  Elastin  is  considered  the 
most  linearly  elastic  biosoild  with  a  Young’s  modulus  of  0.6  MPa  at  100%  stretch  (Park 
and  Lakes  1992).  Elastin  is  a  long,  flexible,  cross-linked  molecule  whose  molecular 
configuration  changes  with  stretch.  Collagen  is  the  basic  structural  element  in  all  tissues, 
with  wavy  fibers  acting  as  mechanical  units  that  dominate  the  force-stretch  response.  The 
initially  compliant  fibers  stretch  and  reorient  under  a  tensile  load  producing  a  nonlinear 
stiffening  response  with  increased  levels  of  stretch  (Fung  1993;  Liu  and  Bilston  2000). 

Holzapfel  suggests  that  soft  biological  tissues  and  solid  polymers  are  unique  in 
that  they  are  the  only  materials  that  experience  a  nonlinear  mechanical  response  to  finite 
strains  relative  to  an  equilibrium  state  (2000).  Since  the  polymer  chain  networks  mimic 
the  assembly  of  the  wavy  cross-linked  collagen,  the  model  used  to  describe  the  behavior 
of  collagen  is  derived  from  models  used  in  rubber  elasticity.  In  particular,  the  3D  network 
model  first  described  by  Arruda  and  Boyce  (1993)  was  used  to  describe  the  collagen  in 
the  parenchyma,  and  a  2D  version  of  this  model,  proposed  by  Prevost,  was  used  to 
represent  the  capsule  (2006). 

The  ground  substance  of  the  ECM,  and  the  cellular  content  of  soft  tissues,  resist 
compressive  forces,  and  maintain  tissue  structure.  Under  a  compressive  force,  these 
molecules  rearrange  to  produce  a  time-dependent  viscoelastic  response.  Therefore,  a 
model  with  a  spring  (to  provide  a  back-stress)  in  parallel  with  a  viscous  component  (to 
account  for  the  time-dependant  viscous  shear  relaxation)  was  used. 

Lastly,  the  extracellular  fluid  in  the  liver  (blood,  bile,  lymph,  etc.)  was  accounted 
for  by  treating  the  liver  like  a  fluid-filled  elastic  sponge.  Using  biphasic  mixture  theory, 
the  fluid  component  responds  to  local  changes  in  volume  by  creating  a  dynamic  pressure 
term  allowing  for  fluid  flow  within  a  porous  elastic  network.  This  flow  in  turn  provides 
resistance  to  changes  in  volume  with  load  over  time.  The  solid  component  was 
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represented  by  an  elastic  spring  whose  bulk  modulus  contributes  to  changes  in  volume  by 
resisting  the  flow  from  the  fluid.  This  component  works  with  the  other  elastic 
components  described  above  to  establish  the  equilibrium  state  of  the  tissue. 

Isochoric  Response  Volumetric  Response 


,ttOr  /Jum>  Gp,  Sp,  m,  fw,  y q  Kg,  KUm,  kq,  M 

Figure  4. 1  Rheological  representation  of  constituve  model  sperated  into  the 
isochoric  (deviatoric)  response  and  volumetric  (hydrostatic)  response.  Material 
parameters  for  the  model  elements  are  listed  below  the  respective  figure. 


4.2  Constitutive  Framework 

The  constitutive  relation  chosen  to  describe  the  liver  is  separated  into  its  isochoric 
and  volumetric  components  and  can  be  represented  by  the  rheological  model  in  Figure 
4.1.  The  total  stress  response  of  the  tissue  ( Ttissue)  under  an  imposed  deformation  (F)  can 
be  obtained  by  summing  the  deviatoric  ( Tisn)  and  hydrostatic  (  Th)  stresses  that  arise  from 
the  isochoric  and  volumetric  component  responses  respectively 

Tissue  =Tiso+Th.  (4.1) 

As  can  be  seen  in  Figure  4.1,  the  total  deformation  gradient  tensor  F  is 
decomposed  into  its  isochoric  component  F  and  its  volumetric  component  Jm  1  (where 
1  is  the  identity  tensor)  and  can  be  defined  as: 


F  =  Jm  IF 


(4.2) 
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where  J  is  the  scalar  volumetric  Jacobian  defined  as  the  current  volume  (V)  divided  by 
the  initial  volume  (Vo),  and  satisfies  the  condition:  det(F)=  1. 

4. 2. 1  Isoch  oric  Network 

Consider  a  multiplicative  decomposition  for  the  isochoric  deformation  where 

F  =  FcFp.  (4.3) 

The  collagen  responds  to  the  Fc  component  of  the  deformation  gradient,  while  the  Fp 
component  is  accommodated  by  the  parenchyma  through  a  dissipative  constitutive 
response.  A  representation  of  the  finite  strain  kinematics  of  this  network  is  shown  in 
Figure  4.2.  Fp  can  be  further  decomposed  into  a  left  stretch  tensor  (Vp)  and  a  rotation 
tensor  (Rp)  so  that 

Fp  =  VpRp  (4.4) 


Configuration  Configuration 

Figure  4.2  Finite  strain  kinematics  for  collagen  network. 


4.2. 1.1  3D  Hyperelastic  Collagenous  Component 

The  formulation  of  the  hyperelastic  collagen  component  of  the  constitutive 
relation  relies  on  the  8-chain  network  model  proposed  by  Boyce  and  Arruda  (2000).  This 
formulation  uses  a  statistical  model  to  represent  the  entropic  state  of  the  individual 
polymer  chains  in  a  network.  A  complete  derivation  of  the  model,  as  well  as  comparisons 
to  Other  such  models,  can  be  found  in  (Arruda  and  Boyce  1993;  Boyce  and  Arruda  2000). 
Although  the  Arruda-Boyce  model  is  derived  from  an  entropic  understanding,  it  has  also 
been  successfully  used  to  represent  the  nonlinear  response  of  biological  networks  where 
the  elastic  response  is  controlled  by  the  internal  energy  (Bischoff,  et  al.  1999).  The  force- 
stretch  relationship  can  be  defined  for  the  individual  fibrils  of  the  collagen  network  by: 
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f^=TL'ii~ )  (4-5) 

b  UJ 

where  A',  is  a  reference  stiffness,  b  is  the  persistence  length  for  the  fibril,  /kand  AL  are  the 
current  stretch  and  locking  stretch  of  the  collagen  fibril,  respectively.  L  is  the  inverse 
Langevin  function  defined  by  (Kuhn  and  Grahn  1942): 

A  A 

—  ,  where  (4.6) 

y 

L(/?)  =  coth/(-l.  (4.7) 

A  simple  representative  cubic  unit  structure  defines  the  geometry  of  the  8-chain 
network  such  that  the  response  of  a  single  chain  (fibril)  can  be  linked  to  the  global 
deformation  of  the  network  due  to  symmetry  (Figure  4.3). 


Bc 


Figure  4.3  The  8-fibril  collagen  network  unit  cell  and  principal  stretches  under 
deformation. 


The  unit  cell  deforms  with  the  principal  stretches,  while  the  junction  point 
remains  centered,  allowing  the  stretch  in  each  fibril  to  be  shared  equally  amongst  all  of 
the  others  via 


In  addition  to  experiencing  the  same  stretch,  each  fibril  also  rotates  towards  the  direction 
of  highest  stretch. 

The  Cauchy  stress  exerted  by  the  collagen  network  is  obtained  in  terms  stretching 
a  collagen  fibril  by  the  left  Cauchy-Green  stretch  tensor  ( Bc ) 

B  =  F  Ft 

c  c  c 


(4.9) 
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Figure  4.4  A  representative  stress-stretch  response  of  the  collagen  netork 
highlighting  a  respone  of  a  single  fibril,  the  initial  and  limiting  shear  moduli,  and  the 
locking  stretch. 


through  the  following  constitutive  law 


T  =  — 
c  J 


f 


Ao 


V 


1.  I 

v  f  J 


M-Ai) 


(4.10) 


where  the  initial  shear  modulus  (po)  and  the  fibril  locking  stretch  (A .£)  are  the  material 
parameters.  The  stress-stretch  response  for  this  relation,  and  a  depiction  of  the  fiber 
orientation  is  shown  in  Figure  4.4.  Since  tissues  cannot  experience  infinite  stress,  and  are 
capable  of  stretches  beyond  the  locking  stretch,  the  slope  of  the  stress-stretch  response 
near  the  locking  stretch  (0.995  X£)  is  used  to  define  a  limiting  shear  modulus  ( pu,n). 

The  collagen  stress  is  considered  to  be  purely  deviatoric  (tr  J^O):  the  collagen  is 
not  allowed  to  account  for  any  of  the  volume  changes  seen  in  the  overall  tissue  response. 
This  “split”  approach  was  taken  due  to  the  observation  that  the  majority  of  the  volume 
change  in  liver  is  due  to  the  movement  of  the  extracellular  fluid,  rather  than  taken  up  by 
the  stretch  of  the  collagen  (Figure  4.5).  An  initial  trial  allowing  the  collagen  to  account 
for  changes  in  volume  proved  problematic  as  the  volumetric  too  dominant.  The  split 
approach  also  has  the  advantage  of  easier  parameter  identification  through  simplified 
initial  calculations. 

4. 2. 1.2  Viscoelastic  Parenchyma  Component 

A  portion  of  the  deviatoric  Cauchy  stress  ( Tc )  derived  above  is  balanced  by  the 
elastic  back  stress  ( Te )  carried  by  the  non-dissipative  component  of  the  parenchyma 
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Figure  4.5  Representation  of  the  assumption  that  the  fluid  component  of  the  model 
takes  up  the  majority  of  the  change  in  volume. 


network.  The  remainder  of  the  stress  (Jv)  drives  the  evolution  of  the  viscous  component 
of  the  deformations  gradient, 

Tc=Te+Tv.  (4.11) 


The  back  stress  element  is  modeled  through  a  linear  elastic  constant  shear 
modulus  ( Gp )  so  that  the  associated  Cauchy  stress  is 


1 

f 

r  trE  "h 

T  =  — 

En  P  1 

e  J 

p 

V 

[p  3  -JJ 

(4.12) 


Here  Ep  is  the  Hencky  strain  associated  with  the  viscous  stretch  (Vp):  Ep=lnVp.  Since  the 
definition  of  a  isochoric  response  requires  that  thedet(F’)  =  det(Fc)  =  det(Fp)  =  1 ,  and 
that  the  det(/?p )  =  det(Kp )  =  1 ,  the  tr(Ep)=Q  and  the  elastic  back  stress  is  reduced  to 

T,=j2GpEp.  (4.13) 

This  response  is  again  purely  deviatoric,  thus  tr(Je)=0.  Te  is  defined  in  the 

instantaneously-unloaded  configuration.  To  properly  prescribe  the  evolution  of  the 
viscous  deformation,  the  back  stress  needs  to  be  pushed  forward  to  the  current 
configuration  by, 

T,C-Fj,FJ.  (4.14) 

The  deviatoric  driving  stress  for  the  viscous  deformation  ( Tv )  is  obtained  from 

TV=TC-Te  (4.15) 

The  magnitude  (z)  and  direction  ( N )  of  the  deviator  of  the  viscous  stress  are  defined  as: 

rp  t  rp  f 

N  =  Tj— 1 )rr  =  JL  ,  where  (4.16) 

PI  V  2r 
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T;  =Tv-itrTvl  (4.17) 

The  viscous  stretching  tensor  ( Dp )  is  constitutively  prescribed  to  be  parallel  to  the 
deviator  of  the  viscous  stress  ( T'v )  through 

Dp=rvN  (4.18) 


where  the  viscous  shear  strain-rate  y v  depends  on  a  power  law  relationship  through  the 


magnitude  of  the  driving  stress  (r)  and  on  material  parameters  Sp  (shear  strength 
modulus)  and  m, 


r  =r0 


\$p  j 


(4.19) 


For  m  =  1 ,  this  relationship  reduces  to  a  linear  viscous  response.  For  m>  1  this  relationship 
can  be  used  to  account  for  the  effect  of  thermally  activated  processes. 

The  rate  of  change  of  the  viscous  component  of  the  deformation  ( F  )  can  be 
expressed  as 

FP  =  LpFp  (4.20) 


where  Lp  is  the  viscous  velocity  gradient  in  the  instantaneously-unloaded  configuration 
and  can  be  expressed  as  the  sum  of  the  stretching  tensor  ( Dp )  and  spin  tensor  ( Wp ), 


Lr=bp+Wr.  (4.21) 

The  stretching  tensor  in  the  instantaneously-unloaded  configuration  is  obtained  by  pulling 
back  the  viscous  stretching  tensor  via 


Dp  =  F~lDpFc .  (4.22) 

The  arbitrary  rotation  associated  with  the  definition  of  the  instantaneously- 
unloaded  configuration  is  eliminated  by  setting  Wp  =0. 

The  stress  from  the  isochoric  network  (7/sy))  is  therefore  defined  as  the  stress  in  the 
collagen  network  (Tc),  which  is  equal  to  the  sum  of  the  stresses  from  the  elastic  ( Te)  and 
viscous  ( Tv )  components  of  the  parenchyma  network  (Tp) 


Chapter  4 


78 


4.2.2  Volumetric  Network 

There  are  two  components  to  the  liver  response  due  to  local  changes  in  volume. 
One  component  considers  an  equilibrium  response  arising  from  the  combined  effect  of 
osmotic  pressure,  volumetric  elastic  response  of  the  solid  matrix,  perfusion  pressure,  and 
resistance  to  extracellular  flow.  The  other  component  arises  from  the  establishment  of  the 
non-equilibrium  pressure  gradient  that  drives  the  fluid  flow  through  the  tissue. 

4.2.2. 1  Equilibrium  Response 

A  simplistic  model  of  the  liver  is  considered  in  which  an  initial  volume  fraction 
(fw)  can  be  associated  with  the  amount  of  free  fluid  in  the  system.  The  remaining  volume 
(1  -fw),  represents  the  solid  phase  and  its  associated  (bound)  hydrated  fluid  tissue 
component.  As  long  as  the  volumetric  strain  (J-l)  is  greater  than  fw,  an  effective 
volumetric  Jacobian  ( J  )  can  be  defined  as 

J *  =  .  (4.24) 

J  w 

The  volumetric  equilibrium  stress  response  (Figure  4.6)  is  then  simply  expressed  as 

<**,=*,  1 4U  (4.25) 

where  Ks  is  a  constitutive  parameter  that  relates  to  the  small  strain  bulk  modulus  (Ko) 
through 


Ks=K0fw  (4.26) 

As  the  volumetric  strain  approaches  the  limit  of fw,  this  formulation  asymptotes  to 
an  infinite  stress  state.  Under  this  condition  it  is  important  to  recognize  that  the  solid 
component  (fs )  is  not  incompressible,  and  can  be  associated  with  a  limiting  bulk  modulus 
Kum-  The  relationship  between  the  hydrostatic  equilibrium  stress  and  the  volumetric 
strain  has  therefore  been  formulated  to  have  an  asymptotic  linear  slope  equal  to  KUm. 
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Figure  4.6  Representation  of  hydrostatic  equilibrium  stress  versus  the  volume 
highlighting  the  equilibrium  (Ks)  and  limiting  (KLim)  bulk  moduli  and  the  amount  of 
free  water  (fw ). 


4.2.2.2  Transient  Response 

Finally,  the  volumetric  transient  response  is  based  on  fluid  flow  through  a  porous 
medium  as  governed  by  a  simple  linear  Darcy’s  Law, 


q  =  -KS/Pfluid 


(4.27) 


where  q  is  the  macroscopic  volume  flow  rate,  k  is  Darcy’s  hydraulic  permeability 
representing  the  resistance  to  flow  under  zero  strain,  and  Pfluid  is  the  dynamic  fluid 
pressure.  To  account  for  variations  in  permeability  with  deformation,  the  permeability  is 
considered  to  be  a  function  of  volumetric  strain  (sv  =7-1)  through  the  relationship  given 
by  (Ateshian,  et  al.  1997) 


K  =  K0 


Ml 

O-tPoks 


(4.28) 


Here  <j)s  and  are  the  respective  solid  and  fluid  contents  of  the  tissue  (where  <j>s+  (j)f= 1), 
Ko  and  </>o  are  the  permeability  and  solid  content  of  the  tissue  in  the  absence  of  volumetric 
strain,  respectively,  and  M  is  the  nonlinear  permeability  coefficient.  k0  and  M  are  material 
parameters  for  finite  deformations. 

The  time  dependent  fluid-pressure  in  the  tissue  can  be  obtained  by  solving  the 
following  equation  at  each  material  point  over  the  liver  domain 
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t/-  P„M=^!Pf,u-C  (4.29) 

_^Lim  J  ^ 

The  volumetric  stress  response  (7/,)  is  therefore  represented  as  the  sum  of  the 
contributions  from  the  hydrostatic  equilibrium  response  {(THeq)  and  the  transient  fluid 
pressure  response  (P fluid), 

Th=crHeql-Pfluidl-  (4-30) 

4.2.3  Hyperelastic  2D  Collagenous  Capsule 

The  liver  capsule  was  modeled  as  a  2D  representation  of  the  3D  hyperelastic 
collagen  network.  This  work  relies  on  the  previous  modeling  work  of  Prevost,  who 
modeled  the  human  chorioamnion  as  an  “idealized”  continuum  membrane  whose  in¬ 
plane  response  is  dominated  by  the  collagen  fiber  network  (Prevost  2006).  A  four-fibril 
network  is  used  as  the  unit  cell,  and  in  this  case  the  collagen  network  is  allowed  to 
respond  to  changes  in  area.  Hence  in  addition  to  the  locking  stretch  ALcap,  initial  (pocap) 
and  limiting  shear  moduli  (/4.imCap),  and  area  expansion  modulus  ( Kcap )  is  also  required  as 
a  material  parameter.  The  Cauchy  membrane  stress  is  given  by, 

Tc  =  JLb^-PU^  +KCv{J-1)Lb  (4.31) 

•'  |_  ^LCap 

where  Bcap  is  the  2D  left  Cauchy-Green  stretch  tensor,  and  fio  is  the  inverse  Langevin 
function  of  (MAuap)-  For  the  capsule,  the  limiting  shear  modulus  was  defined  by  its 
critical  stretch  ( 2unear ).  This  stretch  is  the  point  beyond  which  the  network  starts 
responding  with  a  linear  force-stretch  relationship 

4.3  Physical  Interpretation  of  Framework  in  Relation  to  Observed 
Experimental  Response 

The  constitutive  framework  for  the  response  of  liver  tissue  outlined  in  sections 
4.6.1  and  4.6.2  resulted  in  a  model  for  the  liver  where  the  total  stress  response  from  an 
imposed  deformation  is  a  sum  of  the  contributions  from  the  isochoric  and  volumetric 
responses 

^ tissue  =TC  +  cr Heq  1  —  Pfluid  L  ■  (4.32) 

To  understand  the  physical  meaning  that  correlates  the  biologically  relevant 
material  parameters  to  the  proposed  constitutive  framework,  the  model’s  response  to  a 
step  deformation  is  considered.  The  instantaneous  (t  =  0+),  long  term  (t— >  oo),  and 
transient  (0  <  t  <  oo)  responses  for  different  modes  of  deformation  (change  in  volume, 
pure  shear,  and  indentation)  are  discussed. 

An  instantaneous  change  in  volume  gives  rise  to  very  large  hydrostatic  stresses 
because  Kum  dominates  the  instantaneous  bulk  modulus.  Since  most  biological  tissues 
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have  a  high  water  content  (~70%),  the  modulus  is  expected  to  be  similar  to  those  of 
incompressible  materials,  having  the  same  order  of  magnitude  as  that  of  water  (2.2  GPa). 

Instantaneous  changes  in  shape  (shear  or  isochoric  deformations)  are  accounted 
for  by  an  immediate  shear  response  controlled  by  the  hyperelastic  collagen  network.  The 
collagen  network  provides  the  nonlinear  “stiffening”  character  of  the  tissue  response 
(Figure  4.4).  At  low  levels  of  stretch,  the  curve  has  an  initially  flat  response  described  by 
a  low  initial  shear  modulus  (juo).  The  response  then  sharply  increases  with  stretch  as  the 
limiting  stretched  is  approached  before  becoming  linear  again  with  a  modulus  much 
larger  than  the  initial  shear  modulus  (jium )•  For  a  nonhomogenous  mode  of  deformation, 
like  indentation,  the  material  should  behave  in  a  nearly  incompressible  manner  since  Kum 
is  much  greater  than  the  hyperelastic  shear  moduli  (juo,  pu,n). 

The  dissipative  mechanisms  are  recruited  when  a  constant  level  of  deformation  is 
held  over  time.  These  elements  are  responsible  for  a  portion  of  the  imposed  deformation. 
At  long  times  (t— »oo),  the  material  reaches  an  equilibrium  steady-state  in  which  the 
response  to  changes  in  volume  is  controlled  solely  by  the  equilibrium  bulk  modulus  (Ks). 
The  response  to  change  in  shape  is  accommodated  by  both  the  hyperelastic  collagen 
network  and  the  reconfiguration  of  the  viscoelastic  parenchymal  structure.  This  is 
portrayed  through  a  series  combination  of  the  collagen  moduli  and  a  lower  shear 
parenchyma  modulus  (Gp). 

The  transient  response  of  the  material  is  accounted  for  by  the  material  properties 
of  the  dissipative  components  of  the  model.  The  characteristic  time  for  changes  in 
volume  is  controlled  by  the  permeability  (k)  and  the  length-scale  associated  with  the 
diffusion  paths.  For  a  linear  poroelastic  model  of  a  tissue  specimen  of  height  /  in  confined 
compression,  this  time  can  be  defined  by, 


_  l2 

~  Hk 


(4.33) 


where  H  is  the  confined  compression  modulus.  In  the  case  where  the  elastic  response  of 
the  material  is  nonlinear,  k  is  an  exponential  function  of  volume,  and  the  diffusion  paths 
are  3D  (recall  the  complicated  structure  of  the  liver  vasculature  and  its  hexagonal  lobule 
comprised  of  longitudinal  portal  veins,  hepatic  arteries,  bile  ducts,  and  central  veins  and 
densely  packed  transversely  oriented  fenestrated  sinusoids).  Therefore,  the  relationship  of 
the  volumetric  time  constant  (tvoi)  presented  in  Equation  4.33  is  an  over  simplification. 
However,  it  does  provide  the  correct  dimensional  framework  to  begin  to  interpret  the 
experimental  data.  Note  that  the  size  dependency  of  this  model  feature  cannot  be  replaced 
by  a  simple  viscoelastic  formulation  since  no  length-dependant  behavior  of  the  tissue 
response  is  possible  with  viscoelastic  models. 

Lastly,  the  material’s  transient  isochoric  response  is  described  by  the  nonlinear 
viscous  element  shown  in  Equation  4.19.  This  relationship  is  graphically  depicted  in 
Figure  4.7.  A  linear  response  is  obtained  for  m=  1.  As  the  value  of  m  is  increased,  the 
strain-rate  dependence  diminishes,  and  the  material  is  allowed  to  flow  easily  at  high  rates 
when  the  driving  stress  (r)  exceeds  the  material  strength  (Sp).  These  strain-rate 
independent  behaviors  are  more  indicative  of  biological  materials,  where  it  has  been 
shown  that  hysteresis  is  constant  over  a  wide  variety  of  strain-rates  (Fung  1993),  and  where 
the  tissue  is  filled  with  blood,  a  non-Newtonian  fluid. 
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Figure  4.7  The  effects  of  changing  m  on  the  stress  response  of  the  parenchyma 
viscous  network. 


4. 3. 1  Material  Parameter  Estimation 

The  preceding  sections  sought  to  define  an  appropriate  constitutive  model  and  its 
physiologically  based  parameters  to  appropriately  depict  the  mechanical  response  of  liver 
to  finite  deformations.  Before  implementing  an  iterative  numerical  approximation  scheme 
to  realize  the  values  of  these  parameters,  an  ad  hoc  approach  to  determine  reasonable 
start  values  for  the  parameters  is  presented  based  on  the  literature,  and  features  from  the 
experimental  data. 

•  KLim  is  set  to  0.22  GPa  based  on  the  high  strain-rate  dynamic  testing  reported  by  Saraf 
et  al.  (Saraf,  et  al.  2005) 

•  po  is  bounded  based  on  the  initial  slopes  of  the  ramp  indentation  tests  (1-100  Pa) 

•  pum  is  set  based  on  the  initial  slope  of  the  unloading  curves  during  ramp  indentation 
testing.  The  rapid  change  in  direction  of  the  displacement  is  mainly  accounted  for  by 
the  unloading  of  the  hyperelastic  network.  This  value  was  set  to  1.5  MPa  based  on 
experimental  indentation  data. 

•  XL  is  bounded  based  on  the  curvature  of  the  ramp  indentation  data.  Observing  where 
the  response  turned  nonlinear  approximated  the  locking  stretch  asymptote  (1.01-1.2). 
The  force-displacement  response  of  the  model  is  extremely  sensitive  to  this  parameter 
as  it  drives  the  nonlinear  loading  response,  effects  the  peak  force,  and  dominates  the 
initial  unloading  response  (where  pLim  is  set). 
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•  fw  is  determined  given  the  typical  mass  of  the  liver  (1.1 -1.7  kg)  and  the  amount  and 
densities  of  free  fluid  (namely  blood,  bile,  and  lymph)  stored  within  the  organ  at  any 
time.  Given  values  from  the  literature,  the  fraction  of  free  fluid  in  the  liver  was 
determined  to  be  50-70%  of  it’s  mass,  and  fw  was  set  to  0.5  (Gray  1977). 

•  Gp  was  found  to  be  rather  insensitive  to  ramp  indentation  tests.  Instead  the  stress 
relaxation  data  was  used.  Knowing  that  Kum >  > pum »GP,  this  elastic  shear  modulus 
was  determined  such  that  an  appropriate  steady-state  stress  relaxation  response  was 
achieved  (10-25,000  Pa). 

•  Ks  was  determined  based  on  the  peak  force  of  indentation  ramp  tests,  the  stress 
relaxation  steady-state  response  and  curvature,  and  an  independent  test  conducted 
where  the  mass  of  the  liver  was  measured  over  time  after  a  step  increase  in  perfusion 
pressure.  Prior  to  indentation  testing,  a  harvested  liver  from  an  81  kg  pig  was  allowed 
to  perfuse  for  45  minutes  in  the  ex  vivo  perfusion  system  following  the  protocol  in 
Chapter  3.  The  liver  was  placed  on  a  digital  scale  and  the  I  VC  was  cannulated  so  that 
the  liver  could  drain  elsewhere.  The  hepatic  arterial  perfusion  was  clamped  (to  ease 
the  control  of  the  outflow)  and  the  portal  venous  pressure  was  quickly  raised  76  mm 
inducing  a  step  change  in  pressure  ( dP )  of  0.75  kPa.  The  mass  of  the  liver  was 
recorded  over  time  and  converted  to  volume  ( dV)  (Figure  4.8).  This  procedure  was 
repeated  twice.  The  results  were  fit  to  a  decaying  exponential  whose  constants  ( Veq 
and  Tperf)  were  determined  by  minimizing  the  mean  square  error  and  the  derivative  of 
the  response  to  the  model, 

t 

dV  =  Veq(l -e  l'K"  ).  (4.34) 

Veq  was  found  to  be  3.614x10'  and  5.123  x  10'  m  for  the  two  tests  with  Tperf‘s  of 
32.1  and  44.54  seconds  respectively.  Ks  was  then  estimated  using 

dP 

Ks= - Vg,  (4.35) 

dV  0 

where  Vo  was  the  initial  volume  of  the  liver  prior  to  the  step  change  in  pressure.  The 
results  suggest  an  equilibrium  bulk  modulus  of  37-59  kPa.  Since  the  other  elastic 
moduli  in  the  model  are  also  likely  resisting  the  change  in  volume,  and  recalling  that 
Ks  is  related  to  the  amount  of  free  water  in  the  organ,  Ks  was  bounded  to  be  between  1 
and  50  kPa. 


Chapter  4 


84 


x  10  5 


x  10  5 


Figure  4.8  Change  in  volume  over  time  data  from  filling  test  to  determine  start 
values  for  k() 


k0  was  also  estimated  from  the  results  of  the  swell  test.  Solving  Equation  4.3 1  for  k0, 

l2 


k0  = 


^ petf 


(4.36) 


using  Ks  as  the  moduli,  assuming  a  length-scale  equal  to  \ 4  the  length  of  one  of  the 
liver  lobules  (~0.1  m),  and  using  the  results  of  the  fit  for  Tperf,  initial  estimates  for  k 
were  found  to  be  between  5.8xl0'9  and  6xl0'9  m4/Ns.  Looking  at  the  response  of  the 
model  to  changes  in  k  (namely  peak  force  and  amount  of  hysteresis  in  the  ramp 
indentation  tests)  suggests  that  a: can  range  between  lxl O'6  and  10x1  O'9  m4/Ns. 

•  M  was  determined  from  the  results  of  the  uniaxial  confined  compression  tests  done  on 
cartilage  by  Ateshian  et  al  (1997).  They  showed  that  M  ranges  from  0.4  to  4.3,  and 
that  k  changes  by  10%  when  M  changes  by  >  100%.  Given  the  nonhomogenous 
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indentation  tests,  and  this  insensitivity  of  Mon  the  response,  M  was  set  to  4  assuming 
that  liver  would  fall  on  the  higher  end  of  values  for  cartilage. 

fl  was  set  to  0.01.  Since  the  parameter  for  the  viscoelastic  parenchyma  network  can 

r'n 

be  represented  by  — the  initial  shear  rate  was  set  so  that  flow  would  occur  at  1% 

^  p 

strain/second  when  the  driving  stress  was  equal  to  the  shear  strength  (Sp). 

Sp  and  m  affect  the  strength  of  the  viscous  response  from  the  parenchyma  network. 
This  network  is  assumed  to  be  responsible  for  the  slower  time  constant  of  the 
response.  Therefore,  the  stress  relaxation  data  was  used  to  determine  these  material 
parameters.  Sp  affected  both  the  steady-state  stress  level  as  well  as  the  curvature  of 
the  response,  whereas  m  drastically  affected  the  curvature  of  the  response.  Higher  m ’s 
resulted  in  a  steeper  initial  stress  relaxation.  The  value  of  m  was  set  for  each  liver 
using  the  stress  relaxation  tests,  and  checked  in  the  multiple  rate  ramp  indentation 
tests  to  ensure  that  hysteresis  remained  consistent  with  rate.  Since  the  response  was 
more  sensitive  to  m  than  Sp,  the  range  of  m  was  set  from  1.1-10,  and  Sp  was  more 
loosely  bounded  to  be  between  0.1-100  kPa. 


In  summary,  the  12-parameter  physically-based  constitutive  model  presented  in 
this  section  can  be  reduced  to  five  dependent  material  parameters  (KUm,  pLim,  yvQ ,  fw,  M), 

and  seven  independent  material  parameters  (juo,  Xl,  Ks,  Gp,  Sp,  m,  ko).  The  identification 
of  these  seven  material  parameters  is  done  through  iterative  inverse  finite  modeling  and  is 
discussed  below. 


4.4  Parameter  Identification  via  Iterative  Inverse  FEM 

4.4.1  Numerical  Approximation  of  Constitutive  Model 

Although  indentation  testing  lends  itself  well  for  obtaining  the  mechanical 
response  of  whole  perfused  livers  to  known  forces  and/or  displacements,  the 
interpretation  of  the  results  is  complicated,  particularly  under  finite  deformation.  In 
addition,  the  constitutive  relationship  defined  in  the  preceding  section  (Equation  4.32 
from  Equations  4.10,  4.14,  4.19,  4.25,  and  4.30)  cannot  be  analytically  solved  in  a 
straightforward  manner.  The  two  components  of  stress  giving  rise  to  the  deviatoric 
response  ( Tc ,  Tp)  are  entirely  defined  by  solving  the  initial  boundary  value  problem  by 
following  the  evolution  of  the  deformation  gradient  at  each  material  point.  The 
volumetric  component  can  only  be  obtained  by  solving  the  diffusion  equation  at  each 
material  point  over  the  entire  domain.  The  resulting  nonlinear  multiphysics-based 
problem  cannot  be  solved  without  the  use  of  numerical  approximation  methods  such  as 
finite  element  modeling. 

To  solve  the  inverse  problem  (identifying  the  material  parameters  from  the  known 
force-displacement  response),  the  model  was  implemented  using  commercial  finite 
element  software  (ABAQUS  version  6.5,  Providence,  RI).  The  constitutive  equations 
describing  the  behavior  of  the  2D  liver  capsule  and  3D  tissue  structure  were  coded  as 
user  defined  materials  (UMAT)  in  a  subroutine.  The  implicit  formulation  of  the  time 
integration  used  in  the  subroutine  is  beyond  the  scope  of  this  work,  but  details  can  be 
found  in  Appendix  A  of  Febvay  (2003).  The  implementation  of  the  isochoric  network  is 
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straightforward.  The  diffusion  equation  however  is  realized  to  be  analogous  to  a  heat 
diffusion  equation, 

PcpT  =  kt'V2T  +  rgen ,  (4.37) 

where  p  is  the  density,  cp  is  the  specific  heat,  T  is  temperature,  kt  is  the  thermal 
conductivity,  and  rgen  is  the  heat  source  term.  The  local  equation  for  the  time  evolution  of 
the  fluid  pressure  is  therefore  treated  as  a  coupled  thermal-mechanical  problem. 

4.4.2  Inverse  Finite  Element  Modeling 

Solving  the  inverse  problem  requires  replication  of  the  experimental  boundary 
conditions  in  the  form  of  a  finite  element  model  (Figure  4.9).  Either  the  force  or  the 
displacement  data  are  used  as  prescribed  inputs  to  the  model.  The  constitutive  relation 
and  its  material  parameters  govern  the  model’s  response  to  the  prescribed  loading 
condition.  Initial  estimates  are  given  to  the  material  parameters  and  an  error  estimate 
between  the  model’s  response  and  the  experimental  data  is  determined.  The  parameters 
are  then  updated  in  an  iterative  fashion  until  the  model’s  result  matches  that  of  the  data. 

Given  the  length  scale  of  the  liver  lobe  (-250  x  130  mm)  to  the  size  of  the 
indenter  (6-12  mm  diameter),  the  indentation  tests  were  represented  with  an 
axisymmetric  geometry  assuming  an  infinite  half-space  (Figure  4.10).  The  width  ( w )  and 
height  (h)  of  the  mesh  were  changed  according  to  the  actual  dimensions  of  the  liver.  The 
tissue  was  modeled  as  a  homogenous,  deformable,  solid  with  8-node,  axisymmetric, 
thermally-coupled,  quadrilateral,  biquadratic,  bilinear,  temperature  elements  (CAX8T). 
The  mesh  was  biased  to  have  increased  density  in  the  area  under  the  indenter,  but  kept  as 
coarse  as  possible  to  increase  computation  time.  Comparing  the  results  from  meshes  with 


Figure  4.9  A  schematic  of  the  iterative  algorithm  used  to  solve  the  inverse 
problem  with  a  finite  element  model  (FEM). 
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lOx  finer  resolution  showed  no  noticeable  change  in  reaction  force.  The  total  number  of 
elements  ranged  from  147  to  160  depending  on  the  size  of  the  specimen,  but  each  had  a 
6  mm  x  12  mm  area  directly  under  the  indenter  with  47  elements. 


ho 


Figure  4.10  (Top)  A  schematic  representation  of  the  indentation  experimental 
conditions.  Middle)  The  axisymmetric  meshed  geometry  and  boundary  conditions 
used  in  ABAQUS  to  represent  the  experimental  conditions.  (Bottom)  An  example  of  a 
deformed  mesh  showing  the  reaction  force  on  the  indented  area. 
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The  capsule  was  modeled  as  a  100  pm  thick  homogenous  shell  using  the  skin 
feature  in  ABAQUS.  The  mesh  was  matched  in  number  to  the  tissue  below,  with  3-noded 
quadratic  axisymmetric  thin-shelled  elements  (SAX2).  An  independent  model  was 
created  to  solve  for  the  capsule  parameters  using  the  uniaxial  tension  experiments.  The 
user  defined  material  for  capsule  required  that  an  axisymmetric  condition  be  modeled. 
Therefore  a  tube  in  uniaxial  loading  approximated  a  thin  sheet  in  confined  uniaxial 
tension,  and  a  correction  factor  for  the  force  was  used  to  relate  their  areas 


F 

77  _  model 

cap  ~  0.9 9k 


(4.38) 


Twenty,  2-node,  linear,  axisymmetric,  thin-shelled,  elements  (SAX1)  100  pm  thick  were 
used  to  model  the  capsule.  The  bottom  of  the  capsule  model  was  held  fixed,  the  top  was 
constricted  to  vertical  motions,  and  the  displacement  history  from  the  data  was  used  for 
the  loading  condition. 

Since  the  motorized  indentation  tests  used  a  suction  indenter,  the  complete  liver 
model  was  deformed  by  directly  displacing  the  nodes  equal  to  the  radius  of  the  indenter. 
The  reaction  force  was  taken  from  a  master  node  constrained  to  the  load  nodes.  A  12  mm 
diameter  discrete,  analytical,  rigid,  indenter  with  a  fillet  radius  of  0.4  mm,  was  also 
created  and  the  contact  problem  solved  by  allowing  a  no-slip  condition  (tangential 
friction  set  to  0.9)  under  the  area  of  the  indenter,  and  a  frictionless  condition  on  the 
comers.  Normal  contact  was  defined  using  an  augmented  Lagrange  “Hard”  contact 
pressure-overclosure  for  the  no-slip  area,  and  an  exponential  pressure-overclosure  for  the 
frictionless  section  with  pressure  equal  to  le6  Pa  and  a  clearance  of  4e-5  m.  There  were 
minimal  differences  seen  in  the  force  responses  from  the  model  with  direct  indentation  to 
the  one  using  the  indenter.  Therefore,  the  direct  indention  was  used  since  it  greatly 
decreased  computation  time.  The  6  mm  diameter  indenter  for  the  creep  experiment  was 
modeled,  however,  and  a  frictionless  condition  like  the  one  described  above  was  used. 

The  displacement  data  from  the  motorized  indentation  experiments  were  directly 
implemented  into  the  loading  conditions  for  the  model  using  the  AMPLITUDE  feature  in 
ABAQUS.  The  displacement  history  was  applied  to  the  master  node,  while  the  reaction 
force  was  recorded.  For  the  creep  tests,  an  initial  displacement  was  imposed  on  the 
indenter  to  simulate  the  inertial  loading  of  the  dropped  mass  (-20%  nominal  strain  in 
0.1  sec),  and  then  a  -0.3188  N  load  was  maintained  for  the  duration  of  the  test.  In  all 
cases  the  bottom  of  the  tissue  was  allowed  to  slip  in  the  transverse  direction.  Lastly,  the 
temperature  was  initially  set  to  be  0  everywhere,  and  then  was  inactivated  to  resemble  a 
pressurized  tissue  covered  in  a  non-permeable  membrane. 

4.4.3  Optimization  Technique 

The  seven-dimensional  parameter  space  of  the  constitutive  model  was 
systematically  explored  and  the  results  compared  to  the  data  using  an  optimization 
scheme.  A  link  between  Matlab  (version  7  R14,  Mathworks  Inc,  Natick,  MA)  and 
ABAQUS  environments  carried  out  an  iterative  adjustment  of  the  material  parameters 
according  to  the  Nedler-Mead  simplex  method  (Jordan,  et  al.  2005a;  Lagarias,  et  al.  1998).  The 
objective  function  (0)  to  be  minimized  was  defined  using  the  mean  square  error  between 
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the  data  and  the  model  of  the  force  versus  displacement  (or  time  in  the  case  of  multiple 
indentations)  curve  and  the  time  derivative  of  that  curve 


<D  = 


(F  -  F  ) 

\  data  mod  el  / 


2  A 


+ 


(f  -F  ) 

V  data _ mod  el  ) 


(4.39) 


The  simplex  was  allowed  to  only  sample  the  solution  space  within  the  bounds 
given  for  each  parameter,  and  convergence  was  highly  sensitive  to  the  initial  parameter 
estimation.  Convergence,  or  termination,  conditions  were  set  to  either  a  maximum  of 
1,000  iterations,  or  a  normalized  simplex  diameter  less  than  1  x  104.  If  a  numerical 
solution  was  not  possible,  and  the  model  “crashed”,  an  large  error  (MOO)  was  returned 
and  the  iterations  continued.  An  iterative  solution  for  a  single  indentation  consisting  of 
200  iterations  required  approximately  10  hours  of  computation  time  on  a  3.0  GHz 
Pentium  4  machine  with  2  GB  of  RAM. 


4.5  Results 

4.5.1  Uniaxial  Capsule  Tests 

The  optimization  scheme  was  first  tried  for  the  4-parameter  membrane 
constitutive  model  (Equation  4.31)  on  the  capsule  uniaxial  tension  tests.  Initial  values  for 
the  material  parameters  ( pocap ,  A icaP ,  Aunear,  Kcap)  were  estimated  by  manually  fitting  the 
models  response  to  the  data.  The  optimization  was  then  run  for  each  of  the  9  tests  until  a 
minimum  was  found  (mean  MSE  =  0.00184  in  195-595  iterations).  The  results  for  the 
parameters  are  shown  in  Table  4.1,  and  the  mean  results  are  plotted  against  the  data  for 
each  test  in  Figure  4.11. 

The  mean  and  standard  deviation  of  the  results  suggest  that  the  model  was  more 
sensitive  to  changes  in  the  stretch  parameters  than  the  shear  and  bulk  moduli.  The  model 
response  is  in  excellent  agreement  with  the  data.  The  combined  values  for  the  parameters 
were  used  in  the  indentation  testing  models  below. 


Table  4.1  The  optimized  material  parameters  for  the  liver  capsule  from  uniaxial 
tension  testing 


^-linear 

^-LCap 

I^OCap  [P<0  K 

Cap  [MPa] 

MSE 

Test  1  Mean 
(SD) 

1.023 

(0.021) 

1.040 

(0.02) 

7383.364 

(2830.034) 

1.040 

(0.510) 

0.0022 

(0.002) 

Test  2  Mean 
(SD) 

1.045 

(0.021) 

1.061 

(0.021) 

7385.913 

(2266.511) 

0.671 

(0.249) 

0.00144 

(0.001) 

Combined  Mean 
(SD) 

1.026 

(0.021) 

1.043 

(0.019) 

7081.432 

(2507.71) 

0.948 

(0.483) 

0.00184 

(0.002) 
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Figure  4. 1 1  Results  of  uniaxial  capsule  model  with  mean  parameter  values  and  the 
test  data. 


4.5.2  Liver  Ramp  Indentation  Optimization 

Since  each  parameter  of  the  constitutive  relation  is  sensitive  to  the  type  of  loading 
condition  the  tissue  experiences,  several  deformation  modalities  were  looked  at 
concurrently  to  identity  appropriate  values  for  the  parameters. 

The  data  from  the  multiple  ramp  indentation  and  stress  relaxation  tests  described 
in  Chapter  3  were  used  to  identify  the  seven  parameters  of  the  constitutive  model  for  liver 
defined  in  section  4.3.  Initially,  the  results  of  the  second  indentation  for  each  speed  from 
the  first  block  of  indentations  (rate  order:  2,  0.2,  20,  40  mm/s)  were  used  to  gain  an 
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understanding  of  the  effect  each  parameter  had  on  the  model’s  response.  The  results  of 
the  parameters  were  then  used  and  compared  to  the  stress  relaxation  data.  Fitting 
parameters  to  both  of  these  modalities  in  parallel  provided  confidence  in  the  start  values 
used  for  an  optimization  of  a  full  block  of  multiple  ramp  indentation  tests.  It  was  assumed 
that  the  complete  time  history  of  9-12  ramp  indentations  with  a  120  second  dwell  time  at 
rates  spanning  at  least  2  orders  of  magnitude  would  capture  the  contribution  from  each 
parameter.  However,  since  these  models  were  computationally  expensive  to  run  (>20 
minutes  per  trial),  the  start  values  obtained  from  looking  at  the  stress  relaxation  and 
individual  indentations  at  each  speed  was  necessary. 

The  estimated  start  values  and  allowable  ranges  for  the  parameters  outlined  in 
section  4.3  were  used  to  manually  obtain  good  fits  to  the  stress  relaxation  and  individual 
indentation  data  at  0.2,  2,  and  20  mm/s  for  the  first  liver.  Start  values  were  selected  from 
these  studies,  and  an  optimization  trial  was  run  to  identify  6  of  the  parameters  (m  was  set 
to  2  from  manual  matching  to  the  stress  relaxation  data).  The  force  versus  displacement 
response  for  each  indentation  was  compared  to  the  model’s  response  after  448  iterations 
(Figure  4.12).  The  total  error  was  based  on  the  fit  of  all  three  curves  to  the  data.  The 
evolution  of  the  parameters  and  the  error  are  shown  in  Figure  4.13.  The  resulting 
parameters  were  run  on  a  stress  relaxation  model  and  the  results  compared  to  the  data 
(Figure  4.12).  Due  to  the  poor  fit,  an  optimization  trial  was  then  run  on  the  stress 
relaxation  data  using  the  final  parameters  obtained  for  the  indentation  test  as  start  values. 
This  was  carried  out  for  337  iterations,  and  included  m.  Results  for  this  test  and  the 
evolution  of  the  error  are  shown  in  Figure  4.14.  With  better  confidence  in  the  start  values 
and  with  m  set  to  2.5,  a  final  optimization  was  run  to  identify  the  remaining  6  parameters 
using  the  entire  set  of  9  indentations  and  the  120  sec  dwell  period  for  the  first  block  of  the 
first  liver  comparing  the  force  versus  time  data.  These  models  took  ~20  minutes  to 
complete  per  trial.  Final  parameters  were  obtained  after  64  iterations  (Table  4.2).  The 
resulting  force  versus  time  plots  for  the  block  of  multiple  indentations  across  speed  for 
the  model  and  the  data,  as  well  as  the  error  evolution,  are  shown  in  Figure  4.15. 


Table  4.2  Final  optimized  parameter  values  for  liver  1.  MSE  is  the  mean  squared 
error  between  the  model  response  and  the  data  force  versus  time  curves  for  all  9 
indentations  and  the  120-second  dwell  time.  Units  for  moduli  in  Pa,  k0  is  in  m4/Ns. 


Ks 

gp 

m 

sP 

1^0 

Xl 

K0 

MSE 

25,048 

116.04 

2.5 

13,038 

26.06 

1.039 

6.4 

0.1675 

The  results  show  that  the  model  is  able  to  capture  many  of  the  salient  features  of 
the  time -history  data,  including  an  increased  peak  force  from  the  first  indentation, 
increased  peak  force  with  loading  rate,  decreasing  peak  force  with  multiple  indentations, 
nearly  strain-rate  independent  hysteresis,  and  the  120  second  dwell  time. 

Despite  being  able  to  capture  many  of  the  characteristics  of  the  liver’s  response,  the 
model  fails  to  completely  predict  the  observed  behavior.  There  is  clearly  an 
underestimated  hysteretic  response,  and  an  over  exaggerated  nonlinear  curvature.  The 
predicted  rapid  increase  in  loading  force  during  stress  relaxation  also  disagrees  with  the 
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data.  However,  realizing  that  absolute  parameter  values  will  take  some  time  to  achieve 
(discussed  below),  these  values  are  used  as  estimates  for  the  initial  guesses  of  the  other 
liver’s  material  parameters. 


0.2  mm/s  2  mm/s  20  m m/s 


Figure  4.12  (Top)  Results  showing  the  model  response  (initial  and  final  runs  of  the 
optimization)  to  the  individual  indentation  data  at  three  rates  for  Liver  1 .  (Bottom)  The 
response  using  the  resulting  parameters  compared  to  stress  relaxation  of  the  same  liver 
shown  in  linear  and  log  time. 
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Figure  4.13  The  evolution  of  the  error  and  the  parameter  search  for  the  448 
iterations  comparing  the  model  response  to  the  force-displacement  data  of  individual 
indentations  at  0.2,  2,  20  mm/s  for  Liver  1  (evolution  of  Figure  4.12). 


Chapter  4 


94 


i 


5 

s- 

O 

LU 


°0  50  100  150  200  250  300  350 

Ite  ration 

Figure  4.14  (Top)  The  initial  and  final  optimized  model  responses  compared  to  the 
data  for  the  first  stress  relaxation  test  in  Liver  1  in  both.  (Bottom)  The  evolution  of  the 
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Figure  4.15  (Top)  The  optimized  model  response  to  the  entire  block  of  9 
indentations  for  Liver  1  as  force  versus  time.  (Bottom)  The  evolution  of  the  error. 
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The  third  liver  was  modeled  next,  as  the  second  had  some  inconsistencies  in  the 
data  that  the  model  might  not  be  able  to  capture.  The  data  from  the  second  indentation  of 
each  speed  for  the  first  block  of  12  indentations  was  brought  through  the  same  process 
described  above.  The  first  block  was  chosen  so  the  first  three  speeds  were  done  in  the 
same  order  as  the  first  liver  modeled  (2,  20,  0.2,  and  40  mm/s).  Using  the  values  for  the 
parameters  in  Table  4.2  as  a  guide,  the  model  parameter  values  were  manually  adjusted 
on  the  20  mm/s  indentation  and  the  stress  relaxation  data,  until  a  decent  match  was 
obtained  (Figure  4.16).  The  optimization  across  4  speeds  and  all  7  parameters  was  run  for 
134  iterations.  Although  the  MSE  was  slowly  improving,  the  results  were  converging 
poorly.  Therefore,  the  values  for  the  parameters  at  run  134  were  used  to  fully  optimize 
the  stress  relaxation  response  instead.  The  results  of  this  optimization  (483  iterations)  are 
shown  in  Figure  4.17.  The  single  indentation  at  4  speeds  was  run  again  with  the  newly 
identified  parameters  as  the  start  values  and  m  set  to  the  stress  relaxation  optimized  value 
of  1.18.  This  search  was  stopped  after  158  iterations  as  the  results  were  again  converging 
poorly.  The  force  versus  displacement  data  for  the  4  indentations  is  shown  in  Figure  4.18 
along  with  the  “optimized”  responses  for  each  of  the  trials  run.  The  resulting  values  for 
the  parameters  from  the  stress  relaxation,  and  both  indentation  optimizations  are  shown 
in  Table  4.3. 


Figure  4.16  The  model’s  response  to  the  first  block  of  the  3rd  liver’s  second 
20  mm/s  indentation.  The  red  dotted  line  used  the  parameters  identified  for  Liver  1. 
The  blue  dotted  line  is  a  manual  fit  used  to  obtain  start  values  for  the  optimization. 


Due  to  the  poor  agreement  in  parameter  values  seen  in  Table  4.3,  an  optimization 
for  the  entire  block  was  not  attempted.  However,  a  model  of  the  full  block  was  run  using 
the  mean  values  of  the  parameters,  and  the  results  are  shown  in  Figure  4.19.  Although 
clearly  not  an  exact  fit,  the  model  is  capable  again  of  capturing  the  salient  features  of  the 
data  with  the  exception  of  the  40  mm/s  peak  forces.  Similar  issues  as  were  described  for 
the  first  liver  are  observed.  Namely,  the  model  fails  to  capture  the  true  hysteresis,  and  is 
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unable  to  obtain  an  appropriate  loading  response  for  all  speeds.  In  particular,  the  40  mm/s 
indentations  not  only  had  peak  forces 


Table  4.3  Constitutive  law  parameter  values  for  Liver  3  from  three  different 
optimization  trials.  Units  for  moduli  in  Pa,  k()  is  in  m4/Ns. 
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Trial  2 
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1.18 
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Figure  4.18  Model  responses  for  two  different  optimization  trials  compared  to  the 
second  indentation  of  each  speed  for  liver  3.  Black  is  the  data,  red  is  the  134th  trial  of 
the  optimization  of  all  7  parameters  from  the  manually  fir  start  estimate  in  Figure  4.17 
(MSE  =  0.6742),  and  blue  dashed  is  the  158th  trial  of  an  optimization  with  the 
parameter  results  from  the  stress  relaxation  test  as  start  values  (MSE=1.39). 
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Figure  4. 17  (Top)  The  model  response  to  stress  relaxation  for  liver  3  compared  to 
the  data  with  both  initial  and  optimized  parameters  shown  in  both  linear  and  log  time. 
(Bottom)  The  evolution  of  the  error  function. 
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Figure  4.19  Force  versus  time  data  from  the  first  block  of  data  for  liver  3  and  the 
model  response  using  the  mean  material  parameters  in  Table  4.3. 


Similar  issues  as  were  described  for  the  first  liver  are  observed.  Namely,  the 
model  fails  to  capture  the  true  hysteresis,  and  is  unable  to  obtain  an  appropriate  loading 
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response  for  all  speeds.  In  particular,  the  40  mm/s  indentations  not  only  had  peak  forces 
that  were  20%  of  what  they  should  be,  but  also  were  smaller  than  the  forces  for  the  2  and 
20  mm/s  indentations.  This  suggests  that  perhaps  not  enough  fluid  returned  to  the 
indented  volume  indicating  that  the  permeability  was  too  low. 

4.5.3  Validation 

Although  it  was  already  shown  that  the  results  of  one  liver  might  not  necessarily 
be  true  for  another  liver,  testing  the  predictive  capability  of  the  model  can  still  validate 
the  identified  parameter  ranges.  This  was  done  by  using  the  final  values  of  the 
parameters  in  Tables  4.2  and  4.3  and  modeling  the  response  to  the  creep  response.  Figure 
4.20  shows  the  results  of  this  validation  using  the  data  from  the  first  liver  (recall  that  all 
creep  data  were  similar  in  Chapter  3).  The  material  parameters  for  the  first  liver  capture 
the  initial  dynamic  behavior,  but  do  not  bode  well  over  time.  The  mean  material 
parameters  identified  for  the  third  liver  actually  produce  a  better  fit  to  the  data, 
particularly  in  the  steady-state  case.  This  suggests  that  the  elastic  moduli  are  too  stiff  for 
the  first  liver,  since  they  are  half  the  magnitude  of  those  in  the  third  liver,  and  that  the 
viscous  parameters  for  the  third  liver  are  better  matched  to  the  liver’s  actual  response: 
having  a  softer  shear  stretch  and  lower  exponent  for  the  parenchyma  network. 

Although  not  uniquely  identified,  the  results  of  the  various  optimization  trials 
suggest  a  working  range  for  each  parameter.  These  are  presented  in  Table  4.4. 


Table  4.4  Range  of  material  parameter  values  during  the  optimization  of  the  first 
and  third  livers. 
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Figure  4.20  Model  prediction  of  the  creep  response  of  liver  1  using  the  parameter 
identified  for  liver  1  in  Table  4.2,  and  the  parameters  found  for  liver  3  in  Table  4.3. 


4.6  Discussion 

This  chapter  presents  the  development  of  a  physically-based  nonlinear 
constitutive  model  for  liver,  under  finite  deformations  typical  of  surgical  manipulation. 
Although  others  have  tried  to  model  the  mechanical  response  of  liver  (Brown,  et  al.  2003; 
Carter,  et  al.  2001;  Chen,  et  al.  1996;  Hu  and  Desai  2003;  Kim,  et  al.  2003;  Liu,  et  al.  2004;  Liu  and  Bilston 
2000;  Miller  2000b;  Ottensmeyer  2001;  Picinbono,  et  al.  2001),  no  one  has  captured  the  finite 
deformation,  viscoelastic  response  of  liver  with  a  physically-based  model  like  the  one 
described  here.  Most  models  assume  some  form  of  linearity,  incompressibility,  and  are 
incapable  of  handling  large  strain  indentation  over  time.  The  closet  model  is  one 
presented  by  Snedeker  et  al  for  the  kidney  (2005a).  They  model  the  impact  failure  of  the 
entire  kidney  using  a  multi-constitutive  model  approach.  The  model  represents  the 
capsule,  parenchyma,  calyses,  renal  pelvis,  and  fluid  within  the  kidney.  The  capsule  is 
modeled  with  a  nonlinear  Cowper-Symonds  model  (Snedeker,  et  al.  2005b).  The  response  of 
the  parenchyma  is  modeled  using  a  second  order  Mooney-Rivlin  strain  energy  function 
with  second  order  prony  viscosity.  The  calyses  and  renal  pelvis  are  considered  to  be 
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linear  elastic  fluid  filled  tubes  and  the  blood  vessels  are  treated  as  pressurized  tubes  filled 
with  a  linear  elastic  fluid  (represented  with  a  low  shear  modulus).  Despite  obtaining  good 
agreement  to  their  impact  data,  it  is  unlikely  that  a  model  for  the  kidney  under  such 
conditions  could  represent  the  liver  under  surgically  relevant  deformation. 

This  chapter  described  a  model  and  a  technique  towards  realizing  the  material 
parameters  of  the  liver’s  components.  The  constitutive  framework  developed  shows 
promise  is  representing  solid  perfused  organs  like  the  liver,  under  finite  deformations. 
Shortcomings  of  the  model’s  ability  to  predict  the  true  viscoelastic  response  of  the  results 
obtained  in  Chapter  3,  suggest  that  some  improvements  to  the  model  ought  to  be 
considered. 

One  of  the  strengths  of  this  model  is  the  recognition  that  soft  biological  tissues  are 
compressible.  Although  the  total  stress  can  be  separated  into  components  from  its 
isochoric  (deviatoric)  and  volumetric  (hydrostatic)  responses,  assuming  that  the  changes 
in  volume  come  solely  from  the  fluid  response  is  an  over  simplification.  In  reality,  the 
overall  change  in  volume  should  come  from  an  initial  change  in  volume  from  the  fluid 
moving  out  (the  permeability  network)  and  from  the  collagen  network  losing  water  over 
time.  Allowing  the  collagen  to  contribute  to  the  volumetric  response  would  result  in  an 
additional  material  parameter:  a  prestretch  Ao.  Having  a  pre-stretch  is  relevant  since  the 
fibrous  component  of  biological  materials  are  in  a  state  of  pretension  to  balance  the 
osmotic  pressure  depending  on  the  level  of  hydration  of  the  tissue.  This  prestretch  gives 
rise  to  a  residual  stress  observed  in  other  tissues,  which  is  considered  necessary  for 
successful  modeling  of  finite  deformation  behavior  (Bischoff,  et  al.  2000).  The  curvature  of 
the  force  displacement  response  will  likely  improve  (soften)  by  adding  a  prestretch  to  the 
model,  and  by  allowing  the  collagen  to  change  volume,  albeit  slightly.  Jordan  et  al 
presented  the  effects  of  this  prestretch  in  a  breast  pathology  model  (Jordan,  et  al.  2005a). 

Other  potential  changes  to  the  model  include  adding  a  relaxation  term  to  the 
capsule.  The  creep  response  from  the  model  proved  to  be  much  stiffer  than  that  of  the 
data.  Relaxing  the  parameters  of  the  capsule  helped  achieve  a  more  desirable  response. 
The  results  of  the  uniaxial  tension  tests  to  failure  were  likely  not  enough  information  to 
establish  an  appropriate  model  for  the  capsule.  Stress  relaxation  and/or  creep  tests  on  the 
capsule  should  be  carried  out  to  understand  the  time  dependencies  of  the  tissue  (Prevost 
2006).  The  creep  indentation  test  also  assumed  a  ffictionless  contact  condition.  Dedicated 
tests  in  the  future  should  use  the  motorized  suction  indenter  to  ensure  known  boundary 
conditions. 

The  axisymmetric  assumption  inadvertently  forces  the  fluid  to  stay  contained 
within  the  confines  of  the  given  mesh  geometry.  The  fluid  has  nowhere  to  dissipate  to 
during  long  time-constant  tests,  such  as  the  hour  long  creep  tests.  The  model  therefore 
predicts  higher  steady-state  forces  than  what  is  observed.  Future  models  should 
reconsider  the  geometric  assumptions,  and  better  model  the  flow  path  of  the  vasculature. 

Lastly,  while  the  optimization  results  produced  believable  fits  between  the  model 
and  experimental  data,  the  issues  of  global  convergence  and  uniqueness  of  solution  need 
to  be  addressed.  The  results  presented  here  suggest  that  solving  the  inverse  problem  of  a 
7-dimensional  solution  space  is  not  straightforward,  and  that  the  initial  estimates  for  the 
parameter  values  is  critical  in  ensuring  that  convergence  is  met.  There  are  at  least  two 
ways  in  which  to  address  this  critical  dependence  on  the  initial  parameter  values.  A 
densely  sampled,  exhaustive  search  of  the  7-dimensional  parameter  space,  could  be  used 
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to  understand  the  contour  of  the  solution  space.  Although  extremely  helpful  in  providing 
an  overall  understanding  of  the  parameter  space,  and  in  addressing  the  issue  of 
uniqueness,  the  large  computational  expense  for  performing  such  a  search  is  daunting. 
Another  approach  would  be  to  solve  the  model  analytically  by  simplifying  the  results  of 
the  indentation  tests  with  a  uniaxial  compression  approximation.  A  linear  attempt  is 
presented  in  Appendix  B. 

The  technique  presented  in  this  chapter  also  suggests  that  optimizing  different 
modes  of  time-dependant  deformation  simultaneously  may  be  better  than  doing  them 
sequentially.  To  keep  the  computational  cost  of  this  initial  search  to  a  minimum,  it  is 
suggested  that  a  single  indentation  from  each  speed  of  the  multiple  ramp  tests  and  one 
from  the  stress  relaxation  tests  be  optimized  together  for  all  independent  material 
parameters.  The  error  function  may  also  need  to  be  revised.  Other  measures  for  error 
were  tried  in  an  attempt  to  match  peak  force,  hysteresis,  initial  slope,  and  curvature. 
These  measures  failed  to  improve  the  outcome.  However,  other  measures  of  error  should 
be  considered  that  perhaps  weigh  the  short  and  long  term  response  of  the  stress  relaxation 
more  equally,  etc. 

In  conclusion,  this  chapter  describes  a  physically  based,  physiologically  relevant, 
nonlinear  3D  constitutive  relation  for  the  liver.  The  goal  of  this  thesis  work  was  to  create 
a  “ground-truth”  physically  relevant  model  for  liver.  Simplifications  to  the  model  will  be 
required  for  implementation  into  surgical  simulations  that  require  large  deformations  to 
run  in  real-time.  The  next  chapter  addresses  how  to  evaluate  such  simplifications  by 
assessing  the  speed  versus  accuracy  tradeoff  that  exists  in  real-time  soft  tissue  modeling. 
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Appendix  B 

Linear  Analytical  Approximation  of  the  Constitutive  Model 

In  this  appendix,  the  constitutive  model  presented  in  Chapter  4  is  simplified  to  an 
analytical  linear  approximation.  This  representation  can  be  used  to  ensure  that  the  form  of 
the  constitutive  model  is  correct,  and  to  determine  approximate  start  values  for  the 
material  parameters.  A  nonlinear  analytical  representation  could  obtain  values  for  the 
model’s  actual  material  parameters. 

The  indentation  data  is  represented  with  a  uniaxial  compression  approximation. 
The  volume  of  tissue  under  the  indenter  is  assumed  to  be  the  same  as  a  uniaxial  section 
with  dimensions  equal  to  1.5  times  the  radius  of  the  indenter  ( R ),  and  with  the  same 
thickness  (h)  (Figure  A2.1).  With  this  approximation,  stresses  and  strains  can  be  used 
instead  of  forces  and  displacements,  and  both  sides  of  the  model  can  be  assumed  to 
experience  the  same  strain. 


1.5*R 


Figure  A2. 1  Representing  the  indentation  scenario  (Left)  with 
approximation  (Right). 


a  uniaxial 


The  lumped  parameter  model  is  therefore  simplified  into  the  isochoric  and 
volumetric  parts  with  linear  springs  and  dashpots  replacing  the  nonlinear  and  biphasic 
components  (Figure  A2.2).  The  isochoric  side  can  be  treated  as  a  standard  linear  solid 
viscoelastic  model  with  a  shear  modulus  G,  and  the  volumetric  side  is  a  simple  Voigt 
viscoelastic  model  with  bulk  modulus  B. 
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Isochoric  Response:  G  Volumetric  Response:  B 


Figure  A2.2  Schematic  representation  of  the  simplified  tissue  model. 


—  E0S0  —  CTj  —  (7el  +  <Jvl 

A2.1 

<?el  =Elel 

A2.2 

°vi  =  vA 

A2.3 

£  =  £0  +£l 

A2.4 

Cr2  =  Cre2  +  Crv2 

A2.5 

(7  2  —  E  2^2  +  ^2^2 

A2.6 

£  =  S2 

A2.7 

The  ramp  displacement  waveform  is  approximated  by  a  step  displacement  plus  a 
sinusoidal  waveform  for  ease  of  analysis.  Therefore,  assuming  that  both  sides  of  the 
model  experience  the  same  strain  (Equation  A2.7),  the  strain  and  strain  rate  can  be 
represented  by 

s(t)=  £  -  £ i  sin  cot  A2.8 

o{t )  =  -£tco cos  cot  A2.9 

where  £  is  the  initial  strain,  co  is  the  frequency  of  the  sine  wave  and  t  is  time.  The  stress 
response  for  the  sine  wave  part  is  assumed  to  also  be  a  sine  wave  shifted  in  time  by  phase 
shift  8, 

cj{t)  =  -<7i  sin  (cot  +  5)  A2.10 

&(t)  = -<7iCocos(cot  +  8)  A2.ll 

while  the  step  response  is  assumed  to  have  an  exponentially  decaying  stress  response. 
Substituting  equations  A2.8-1 1  into  equations  A2.1-6  results  in  the  equation  of  motion 
for  each  side 
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a"(,)=s‘  \jffk  1_* 


A2.12 


(E,+Ex)2  +  co1  ri\ 


[sin  cot{Ex  (E0  +  Ex )+  co2p\  )+  cos  cot{E  Opxco^[ 


<j2 (t)  =  -E2si  - ei [sin cot(E2)+ cos a>t{r/2co )] 


A2.13 


The  right  hand  term  of  each  equation  of  motion  is  the  response  from  the 
sinusoidal  input,  and  these  can  be  used  to  determine  the  shear  and  bulk  storage  (in  phase) 

and  loss  (out  of  phase)  moduli  (Gs,  Go,  Bs,  BD)  and  the  shear  and  bulk  complex  moduli 

*  * 

0 G,B ) 


Eq(ei(E0  +  £) )+  a)2 i{] ) 
(e0  +Ej y  +m2f]2I 


A2.14 


Gd  = 


(e„ +Ej y +co2f]'I 


A2.15 


G  =  iGS  +  iGD  ) 


A2.16 


E s  —  E2 


A2.17 


Bd  =  cop, 


A2.18 


B  -  (Bs  +  iBD  ) 


A2.19 


Assuming  uniaxial  compression,  an  overall  elastic  moduli  ( E  )  for  the  model  can 
be  represented  by 


E *  =  98  G  =  £  +iED 
G  +3B  s 


A2.20 


where  the  storage  (Es)  and  loss  (Ed)  moduli  for  the  whole  model  can  be  determined  by 
substituting  equations  A2.14-A2.19  into  A2.20. 


9  E0(num) 


A2.21 


'j2  +  orp2  \E;  +  a>2p2 


)+E0(E2(e2  +3ElE2  +co2p2)+3co2Elp2))  A2.22 
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den  = 


A2.23 


9 (e2  +  a2ril  \e22  +  co2ri22 )+  6E0 (e2  ( E2  +  3ElE2  +  co2jj2  )+  3co2E1t]2  )+  E2 ((£,  +  3 E2  )2  +  to 2  (p,  +  3tj2  f ) 


E  =  +3 tj2)))  A2  24 

D  den 

At  high  and  low  frequencies  the  loss  modulus  tends  towards  zero  while  the  storage 
modulus  reduces  to 


lim  (ES)^>3E0  A2.25 

(D-> 00  V 


lim  fe)-> 


9E0E1E2 

3E  jE  2  +E0(E1  +  3E  2) 


A2.26 


The  tangent  of  the  phase  shift  is  taken  to  be  the  loss  modulus  divided  by  the  storage 
modulus 


t  c=  _  coEjQE\r^  +  ?/2  (eI+cq\  fa  +  3rj2 )))  _ 

3(^!2  +®2712X^2  +®2/722  )+£0(^2(£12  +3£t£2  +®27/12)+3  «2^722) 

Assigning  values  for  the  spring  and  dashpot  constants,  the  storage  modulus,  loss 
modulus,  and  tan 5  (Equations  A2.21,  A2.24,  A.27)  are  plotted  versus  frequency  in  Figure 
A2.3.  The  values  for  the  constants  are:  E0  =  5,131.4  Pa  (from  a  spring  constant  Ko  =  50 
N/m),  Ei=  513.15  Pa  (from  a  spring  constant  Ki  =  5  N/m),  E2=  2,565.7  Pa  (from  a  spring 
constant  K2  =  25  N/m),  q  1=  10,263  Pa  s  (from  a  dashpot  constant  bi  =  100  Ns/m),  and 
r|2=  1,026.3  Pa  s  (from  a  dashpot  constant  bi  =  10  Ns/m). 

Given  Figure  A2.3  and  the  equations  presented  above,  the  constants  for  the  model 
can  be  determined.  A  more  robust  interpretation  of  the  model  would  provide  a  better 
estimation  of  the  model  parameters  by  comparing  the  predicted  moduli  to  the  moduli 
from  the  single  ramp  rate  data  collected  in  Chapter  3.  Changes  to  this  analytical  model 
would  include  using  the  appropriate  nonlinear  springs  and  dashpots  and  using  the 
combined  modulus  for  indentation  in  place  of  Equation  A2.20  (although  meant  for 
infinitesimal  strains) 


*_  12G*(-3B*  +G*)2 
~(95*-4G*)(55*+G*) 


A2.28 
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Figure  A2.3  The  storage  modulus  (Top),  loss  modulus  (Middle),  and  tan8  (Bottom) 
of  the  linear  model  across  frequency.  The  different  curves  represent  different  values 
for  the  isochoric  dashpot  [units  for  b=  Ns/m]. 
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ABSTRACT 


The  ability  to  accurately  model  mechanical  properties  of  soft  tissues  plays  a  fundamental  role  in  the  future 
advances  in  surgical  simulation,  planning,  and  assistance.  While  a  physically  based  mathematical  model 
capable  of  capturing  the  complex  behavior  of  soft  tissues  has  been  formulated,  the  identification  of  unique 
model  parameters  remains  a  challenge.  The  goal  of  this  project  is  to  combine  traditional  indentation  testing 
with  three-dimensional  ultrasonic  imaging  to  capture  the  mechanical  response  of  in  vivo  soft  tissue  and,  with 
this  information,  uniquely  identify  a  set  of  optimal  model  parameters.  This  technique  is  currently  being 
developed  for  modeling  the  response  of  liver  tissue  but  will  be  extended  in  the  future  to  other  soft  tissues, 
including  those  of  spleen,  kidney,  breast,  and  brain.  The  development  of  this  technique  is  a  key  step  toward 
image-driven  tissue  modeling  and  characterization. 


METHODS 


Experimental  Approach 

Characteristic  force-displacement  relationships  are  obtained  from  an  indentation  of  porcine  liver  while  imaging  with  3D 
ultrasound  (Philips  SONOS  7500  Live  3D  Echo,  Philips  Medical  Systems,  Andover,  MA,  USA),  providing  the  ability  to 
estimate  the  complete  three-dimensional  deformation  field  of  the  sample  under  indentation. 


The  stress-strain  characteristics  of  biological  tissues  are  known  to  be  largely  nonlinear,  rate  dependent,  and  dissipative. 
Our  soft  tissue  model  [2,3]  consists  of  a  collagen-elastin  structural  component  in  the  form  of  an  extended  Arruda-Boyce 
hyperelastic  law,  a  cellular  interaction  modeled  as  a  nonlinear  viscous  solid,  and  an  organ  vascular  perfusion  modeled  by 
a  pressure-driven  diffusion  process.  The  response  of  the  model  is  fitted  to  in-vitro  porcine  liver  indentation  experiments 
through  iterative  finite-element  modeling. 

•  iterative  FEM  parameter  estimation  framework  is  used 
to  match  the  experimental  response  to  a  mechanical 
organ  model 


•  material  parameters  are  found  by  minimizing: 

O  (A0,XL,jU,k,  K)  =^Yf(Fmodel(di’A0,AL,M,k,  $-Fexperimtnt(dl')f 


I 

-  4.=.  ^  _ 

d 

1 

-1  Etf  1 _ .  -j* 

3  T  (MSE) 

Jd, i  -m 

Optical  Flow  Formulation 


-  - 


t  t  +  st 

The  goal  of  optical  flow  tracking  is  to  estimate  displacement  of 
every  voxel  between  data  frames  collected  at  times  t  and  t  +  St. 


•  Under  intensity  (£)  constancy  assumption,  the  displacement 
(u,v,  w)  of  a  voxel  from  time  t  to  t+St : 

I(x  +  uSt,y  +  vSt,z  +  wSt,  t  +  dt)-  I(x,y,z,  t) 


•  Resulting,  upon  Taylor  series  linearization,  in  the  optical 
flow  constraint : 


dt  dt  dt  dl  ^ 

- u- 1 - V  H - W - =  0 

dx  dy  dz  dt 


•  n  voxels  results  in  n  equations  and  3*n  unknows  and 
requires  additional  constraint 


One  of  the  limitations  of  indentation  testing  methods  is  the  area  of  uniqueness  and  sensitivity  of  model  parameters  to 
experimental  observations.  In  order  to  improve  the  tissue  testing  process,  we  present  a  3D  ultrasound  imaging  technique, 
which  captures  the  experimental  tissue  deformation  during  indentation.  We  employ  computer  vision  techniques  to  track 
tissue  deformation  with  a  novel  optical  flow  algorithm.  Standard  optical  flow  algorithms  are  not  explicitly  consistent  with 
solid  organ  mechanics  and  tend  to  bias  their  deformation  estimate  toward  the  motion  of  an  idealized  incompressible  linear 
solid.  Our  work  has  addressed  this  limitation  by  a  development  of  a  regularization  framework  capable  of  coupling  optical 
flow  equations  to  fundamental  soft  tissue  biomechanics. 


Horn  and  Schunck  Approach  [4] 

Find  voxel  displacements  (u,v,w)  that  minimize: 


'V(cc)  =  jjj  (yV{)(x,y,z,t)  +  cc'V s(x,y,z,t))dxdydz 


¥„«( E,u+E,v+E,w+E,Y 

Resulting  Laplacian  smoothing  operator: 


V2«  =  0 


FEM-Constrained  Optical  Flow 

Optical  flow  solved  as  a  mechanical  system: 


Per-node  least-square 
estimation  of  (u,v,w) 
and  (cx,  cy,  cj 


y  ~\ 

Full  displacement  field  from 
mechanical  FEM  model 


Smoothing  operator  for  linear  elastic  materials: 


(l+v)(l-2v) 


V(V-w)+ 


2(1  +  v) 


-V  ti  =  0 


Mechanics  /  Optical  Flow  Coupling 


Solving  the  Optical  Flow  Field 


■  linear  springs  are  attached  to  the  mechanical 
model,  representing  local  optical  flow  computed  via 
modified  Lucas-Kanade  method  [5] 

■  the  global  optical  flow  is  solved  as  a  standard 
mechanics  problem  with  a  commercial  FEM 
package 

■  resulting  optical  flow  field  reflects  proper  body 
mechanics  for  a  chosen  material  law 


Local  motion  estimates  (spring  parameters): 


Nnoie{x„y„zt) 

OX 

d/(x,.,y,.,z.) 

^node(X„y„Zt) 

dy 

dz 
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'  „  ,  sdl{xt,ytiz,)  ' 

dt 

M 
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M 

M 

dx 

dy 

dz 

-w- 

Nmde{xn,yn,z n)dI^y-^ 
mdA  "  y"  dt 

Where  Nnode(x,y, z)  are  linear  element  shape  functions  node  ^  Ar  ^  j5/(xI.,y,.,zi)| 

and  (cx,Cy,cJ  are  confidence  levels  computed  as:  Cx  —  | 


FEM  implementation 

■  5-parameter  poro-elastic  constitutive  law  [2,3] 

■  Linear  tetrahedral  mesh 

■  Fluid  flow  modeled  with  temperature- 
displacement  coupling 

■  Finite-deformation  formulation 

■  Frictionless  indenter-tissue  boundary  condition 

■  Modeled  in  ABAQUS  6.5-1  (HKS,  Providence,  Rl) 


RESULTS  &  FUTURE  WORK 


Generation  of  Synthetic  Data 


Synthetic  Data  Generation 


■  serves  as  "gold  standard"  for  motion 
tracking  validation 

■  3D  block  of  ultrasound  texture  from 
liver  parenchyma 

■  texture  block  is  deformed  by  uniaxial 
compression,  assuming  liver-like  organ 
mechanics 


Motion  Tracking  Performance 

■  mechanical  regularization  provides 
superior  accuracy  to  Horn  &  Schunck 

■  mechanical  regularization  has  superior 
convergence  characteristics 

■  both  techniques  perform  well  under 

noise  (10%  max-intensity  white 

noise) 


Future  Work 


■  application  of  the  motion  tracking  algorithm  to  soft  tissue  characterization 


■  investigation  of  parameter  sensitivity  under  various  imaging  and  testing 
conditions 


ACKNOWLEDGMENTS 


Funding  by:  U.S.  Army  Medical  Research  and  Material  Command,  contract  number  DAMD  17-01-1-0677. 

References: 

[1]  Kerdok  AE,  Jordan  P,  Liu  Y,  Wellman  PS,  Socrate  S,  Howe  RD.  Identification  of  Nonlinear  Constitutive  Law 
Parameters  of  Breast  Tissue.  Proc.  of  the  2005  Summer  Bioengineering  Conference,  ASME,  2005. 

[2]  Socrate  S,  Boyce  MC.  A  constitutive  model  for  the  large  strain  behavior  of  cartilage.  In  Proceedings  of  2001 
Bioengineering  Conference,  volume  50,  pages  597-598.  ASME,  2001. 

[3]  Febvay  S,  Socrate  S.  Biomechanical  modeling  of  cervical  tissue:  a  quantitative  investigation  of  cervical  tunneling.  In 
ASME  Int.  Mechanical  Engineering  Congress  and  Exposition  (IMECE),  2003. 

[4]  Horn  BPK,  Schunck  BG.  Determining  optical  flow.  Art.  Intelligence,  16(1 -3):  186-203,  August  1981. 

[5]  Lucas  BD,  Kanade  T.  An  iterative  image  registration  technique  with  an  application  to  stereo  vision.  Proc.  of  the  7th 
Inter.  Joint  Conf.  on  Artificial  Intelligence  (IJCAI  ’81),  pages  674-679,  April  1981. 


